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INTRODUCTION 


A.  STATEMENT  OF  THE  PROBLEM 

The  Army  has  a  need  for  information  about  the  distribution  of  smokes  and  chemical  agents 
in  the  atmosphere,  and  their  changes  with  time.  With  the  increasing  use  of  laser-guided 
weaponry  and  detection  systems,  there  is  a  great  need  for  quantitative  measures  of  the  probability 
that  a  clear  line-of-sight  will  be  found  through  an  obscurant  cloud.  Chemical  and  biological 
defense  planning  also  needs  information  about  the  likelihood  that  critical  concentrations  will  be 
exceeded.  This  is  the  problem  addressed  by  the  research  reported  here. 

FYactical  treatments  of  the  smokes  and  gasses  in  the  atmosphere  have  largely  been  limited  to 
dealing  with  “average”  or  “most  probable”  results,  i.e.,  smooth  distributions.  Any  observation  of 
atmospheric  plumes  suggests  that  such  a  picture  does  not  accurately  portray  conditions  at  a  single 
moment  in  time.  Considerable  structure  and  numerous  inhomogeneities  are  usual  features  of 
such  instantaneous  plumes.  This  is  not  important  for  some  applications,  but  in  many  instances 
it  is,  (e.g.,  when  we  wish  to  know  the  likelihood  of  being  able  to  see  through  a  plume,  or  of 
encountering  short-term  concentrations  above  some  critical  threshold). 

There  is  reason  to  believe  that  the  inhomogeneities  in  smoke  plumes  are  “fractals.”  Some  of 
these  reasons  are  reviewed  later.  In  the  work  reported  here,  I  have  identified  suitable  fractal 
analysis  techniques  and  observations,  and  adapted  and  applied  those  techniques  to  the  observa¬ 
tions.  Some  of  the  techniques  and  results  suggest  that,  in  the  future,  they  will  lead  to  a  much 
better  understanding  of  atmospheric  turbulent  processes. 

B.  EVOLUTION  OF  THE  RESEARCH  APPROACH 

Before  proceeding  to  review  why  fractal  concepts  are  important  to  the  description  of  the 
inhomogeneous  distributions  of  kinematic  properties  and  the  mixing  of  materials,  I  will  briefly 
describe  the  evolution  of  the  approach  that  was  taken. 

This  work  started  with  the  intention  of  developing  and  applying  methods  for  calculating  the 
fractal  dimension  of  atmospheric  fields.  As  the  project  proceeded,  it  became  apparent  that  fractal 
dimension  provided  only  limited  information  and  that  the  techniques  used  to  estimate  fractal 
dimension  generated  much  valuable  information  that  was  discarded.  For  example; 

•  Simple  box -counting  methods  do  not  use  the  spatial  relationships  between 
boxes  containing  members  of  the  set  at  different  scales. 

*  The  Fourier  methods  do  not  use  the  distributions  of  amplitudes  in  the  Fourier 
plane  to  infer  anything  about  changes  in  anisotropy  with  scale. 

Other  factors  were  also  troublesome.  First,  all  the  methods  dealt  only  with  the  distribution 
of  scalars  in  space,  or  with  the  distribution  of  members  of  some  specified  set  (e.g.,  scalars  with 


I 


values  in  a  specified  range).  Another  problem  area  was  that  fractal  dimension  was  of  rather  lim¬ 
ited  use.  Finally,  there  didn’t  seem  to  be  much  connection  to  the  physics  of  the  processes  being 
studied,  although  there  appeared  to  be  some  good  reasons  to  think  that  fractals  would  be  inherent 
to  many  atmospheric  processes  (e.g.,  the  space  fillingness  of  turbulence  and  the  cascade  of 
energy  to  ever  smaller  vortices). 

At  about  this  point  J.G.  Jones  sent  a  preprint  of  the  paper,  “Multi-Resolution  Analysis  of 
Remotely-Sensed  Data.”  This  paper  (eventually  published^ — Jones  et  al,  1991)  described  how  to 
use  defined  structures  in  a  two-dimensional  scalar  field  as  the  sets  from  which  fractal  dimension 
was  calculated.  Although  still  limited  to  scalars  and  two  dimensions,  it  was  clear  that  the  tech¬ 
nique  could  readily  be  extended  to  three  dimensions.  Furthermore,  it  soon  became  obvious  that 
extension  of  the  approach  to  two-  and  three-dimensional  vector  fields  was  also  possible.  This 
multiresolution  feature  analysis  and  other  techniques  are  described  in  Section  11. 

With  this  new  approach,  it  was  theoretically  possible  to  select  “physically  significant  fea¬ 
tures,”  and  use  them  in  the  analyses  to  learn  more  about  the  processes.  It  also  appears  that  use  of 
such  features  might  significantly  improve  the  realism  of  simulated  inhomogeneous  distributions. 
Of  course,  we  must  be  able  to  define  “significant  features,”  the  obvious  line  of  reasoning  is  to 
devise  a  method  whereby  the  data  themselves  would  determine  what  was  important.  Statisticians 
frequently  use  principal  component  analysis  for  this  purpose,  and  Lorenz  (1956)  had  successfully 
used  a  closely  related  approach  to  define  patterns  of  pressure  and  temperature  variability  for  sta¬ 
tistical  weather  forecasting.  Lorenz  called  the  patterns  that  account  for  the  most  variance  in  the 
data  Empirical  Orthogonal  Functions  (EOFs). 

The  EOF  technique  has  the  following  advantages: 

•  The  main  patterns  of  variability  are  defined  by  the  data  themselves,  thereby 
giving  reason  to  believe  that  they  are  “physically  significant.” 

•  Because  the  patterns  are  linearly  independent  (orthogonal),  the  original  pat¬ 
terns  of  variability  (and  their  spatial  derivatives)  can  all  be  approximated  by 
linear  combinations  of  features. 

It  also  turns  out  that  EOFs  may  reveal  artifacts  in  the  data  introduced  by  the  instrument  used  to 
collect  the  data. 

In  addition  to  the  Jones  et  al.  (1991)  analysis  technique  that  was  reported  during  the  course 
of  the  project,  new  types  of  data  also  influenced  the  research  direction.  Two  data  sets  were  par¬ 
ticularly  important.  First,  was  a  set  of  detailed,  three-dimensional  wind  observations  provided  by 
Schneider  (1991 ),  and  second  were  detailed  estimates  of  transmittance  through  a  smoke  plume 
presented  as  two-dimensional  imagery  (Bleiweiss  et  al.,  1991).  Transmittance  inhomogeneities 
through  a  plume  represent  an  imponant  class  of  problem  for  the  Army.  Unfortunately,  the  data 
were  obtained  only  a  few  months  before  the  conclusion  of  the  project,  preventing  as  comprehen¬ 
sive  an  analysis  as  desired. 

Considerable  effort  was  spent  in  developing  vector-analysis  techniques  for  application  to 
the  detailed  wind  observations  in  the  belief  that  an  understanding  of  the  atmospheric  motions 
would  provide  a  basis  for  describing  the  inhomogeneities  of  the  scalar  distributions.  While  this 
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still  is  a  reasonable  expectation,  it  became  apparent  while  developing  the  vector-analysis  tech¬ 
niques  that  this  approach  would  require  more  effort  than  could  be  expended  as  part  of  the  cuirent 
project;  therefore,  the  focus  returned  to  the  analysis  of  scalar  distributions.  Although  the  main 
focus  of  this  report  is  on  the  distribution  of  scalars  in  the  atmosphere,  some  of  the  vector-analysis 
approaches  and  preliminary  results  are  also  discussed. 

Some  conclusions  can  be  drawn  from  the  studies  reported  here,  but  to  a  large  extent  they  are 
only  preliminary.  Some  of  the  techniques  appear  promising,  but  have  not  been  refined  and 
applied  rigorously  enough  at  this  stage  to  realize  that  promise.  Thus,  some  of  the  conclusions 
reported  at  the  end  of  this  report  take  the  form  of  suggestions  for  refinements  and  future  applica¬ 
tions  of  the  techniques  described  here. 
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n  REVIEW  OF  FRACTAL  CONCEPTS 
AND  THEIR  APPLICATIONS  TO  THE  ATMOSPHERE 


This  section  is  in  large  part  derived  from  a  review  published  earlier  (Ludwig,  1989),  but 
inasmuch  as  this  is  a  rapidly  evolving  topic  it  has  been  necessary  to  include  some  more  recent 
material.  The  purpose  of  the  review  is  to  provide  definitions  for  fractal  dimension  and  descrip¬ 
tions  of  various  methods  by  which  it  can  he  calculated.  Another  important  goal  is  to  provide  a 
basis  for  understanding  why  the  concept  of  fractals  has  relevance  to  the  understanding  of  atmo¬ 
spheric  inhomogeneities. 

The  word  “fractal”  only  appeared  about  16  years  ago  (Mandelbrot,  1975).  The  widespread 
use  of  the  word  in  popular  and  technical  writing  in  recent  years  suggests  that  either  fractal  con¬ 
cepts  have  general  applicability,  or  that  fractals  have  taken  on  something  of  a  cult  status.  The 
fact  is  that  both  reasons  apply:  Concepts  related  to  fractals  do  provide  good  mathematical 
analogs  in  fields  ranging  from  cosmology  to  geography  to  population  biology  and  fluid  mechan¬ 
ics  (Mandelbrot,  1983);  unfortunately,  there  may  also  have  a  been  a  tendency  to  think  that  this 
widespread  applicability  meant  that  a  “universal  truth”  had  been  found.  Atmospheric  science  is 
among  those  places  where  the  fractal  has  found  a  home.  The  ability  of  fractal  concepts  to 
describe  a  wide  variety  of  conditions  and  processes  is  tantalizing,  but  to  date,  the  practical,  pre¬ 
dictive  results  have  been  few. 

This  review  is  intended  to  provide  a  qualitative  understanding  of  important  concepts  and  an 
introduction  to  the  simplified  version  of  the  mathematics  underlying  these  concepts.  Examples 
of  how  the  concepts  apply  to  the  atmosphere,  and  how  some  relate  to  classical  approaches  are 
also  included. 

A.  GENERAL  CONCEPTS 

One  of  the  more  comprehensible  pans  of  many  discussions  of  turbulence  is  the  following 
piece  of  doggerel  credited  to  Richardson  (1922): 

Big  whorls  have  little  whorls, 

Which  feed  on  their  velocity; 

And  little  whorls  have  lesser  whorls. 

And  so  on  to  viscosity. . . 

Those  four  lines  provide  a  vivid  picture  of  turbulence;  they  also  provide  a  good  staning 
point  for  visualizing  fractals.  If  we  examine  the  streamlines  associated  with  Richardson’s  big  and 
little  whorls,  we  find  that  those  streamlines  have  their  corresponding  large,  little,  lesser,  and  least 
w}'3gles — down  to  where  “viscosity”  stops  the  process.  We  can  magnify  and  find  wiggles  super¬ 
imposed  on  wiggles  already  there.  This  is  one  type  of  fractal.  In  this  case,  we  have  a  collection 
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of  points  that  fonn  a  line  with  wiggles  and  roughness  over  a  large  range  of  scales.  The  set  of 
points  could  also  be  part  of  a  surface  with  roughness  elements  at  every  scale,  or  a  collection  of 
separated  points  that  are  arranged  in  clusters  of  clusters  of  clusters.. . 

Many  lines  and  surfaces  in  nature  are  fractals.  Coastlines  have  the  property  of  being 
“rough”  over  a  wide  range  of  .scales  (Mandelbrot,  1983;  Feder,  1988).  Lovejoy  (1982)  found  thai 
cloud  and  rain  areas  have  these  properties  over  more  than  four  orders  of  magnitude;  Rys  and 
Waldvogel  (1986)  looked  at  hailstorm  perimeters  and  obtained  similar  results  down  to  sizes  of  a 
few  kilometers,  but  found  smoother  outlines  at  smaller  scales. 

Fractals  are  rough  geometric  shapes:  their  roughness  is  qualitatively  similar  at  all  scales; 
the  mea.sure  of  roughness  is  fractal  dimension.  Later  1  discuss  several  quantitative  mathematical 
expre.ssions  that  can  be  used  to  determine  fractal  dimension;  here  I  limit  the  discussion  to  some 
easily  understood  subjective  attributes  associated  with  fractal  dimension,  such  as  the  already 
noted  relation  between  fractal  dimension  and  perceptions  of  roughness.  Pentland  (1984)  reported 
that  when  people  ranked  images  of  surfaces  in  order  of  their  perceived  roughness,  the  rankings 
corresponded  to  fractal  dimension,  with  the  roughest  surfaces  having  the  highest  fractal 
dimension. 

Space-fillingness  is  another  useful  subjective  concept  associated  with  fractal  dimension.  A 
smooth  line  is  confined  to  one  dimension,  a  smooth  surface  to  two.  A  fractal,  with  its  wiggles  on 
wiggles  on  wiggles,  begins  to  infringe  on  other  dimensions.  Topologically,  the  fractal  line  or 
fractal  surface  may  have  only  one  or  two  dimensions,  but  the  wiggles  make  them  substantively 
different  from  the  smooth  shapes.  Fractal  dimension  provides  a  quantitative  measure  of  the 
degree  to  which  these  structures  fill  the  physical  space  beyond  their  topological  dimension. 

The  concepts  discussed  above  can  be  extended  to  any  number  of  dimensions,  but  consider¬ 
able  care  must  be  taken  when  the  dimensions  represent  different  physical  properties  (e.g.,  a  mix 
of  pressure,  density,  and  wind  components).  Even  in  two  dimensions  (e.g.,  a  time  series  of 
aerosol  concentrations  at  a  point),  Mandelbrot  (1985)  warns  that  there  can  be  discontinuities  in 
certain  measures  related  to  fractal  dimension. 

B.  MATilEMATIC.AL  RELATIONSHIPS 

Several  different  definitions  of  fractal  dimension  are  found  in  the  literature.  They  illustrate 
different  physical  properties  or  processes  that  are  related  to  fractal  dimension.  The  first  two 
methods  for  defining  and  calculating  fractal  dimansion  discussed  below  follow  from  the  fact  that 
fractal  dimension  is  a  measure  of  space-filling  properties.  The  other  definitions  are  related  to  the 
fact  that  fractals  are  associated  with  fields  that  are  self-similar  (i.e.,  scaling)  over  a  wide  range  of 
scales. 

1.  Box  (or  Cell)  Counting 

Mandelbrot  (1975)  defined  fractal  dimension  through  its  relation  to  the  space-filling  con¬ 
cept.  He  assumed  a  continuous  distribution  of  some  parameter  so  that  there  are  isosurfaces 
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containing  all  the  points  with  the  same  parameter  value.  He  then  invoked  a  concept  commonly 
used  for  measuring  metric  properties  by  counting  the  number  of  elementary  units  (e.g.,  small 
units  of  area  or  volume)  required  to  cover  all  the  points  in  the  set  Restated,  a  metric  property, 

M,  is  estimated  by  counting  the  number  (N)  of  measurement  units  of  linear  size  X  that  are 
required  to  cover  the  set  of  points  being  measured.  For  simple  shapes,  the  measure  is  related  to 
the  dimension,  D,  of  the  shape  and  the  size  of  the  measurement  unit: 

M  =  NXD  .  (1) 

Mandelbrot  (1983)  has  noted  that  the  usual  smooth  shapes  of  Euclidian  geometry  result  in  inte¬ 
ger  values  for  D,  but  many  natural  shapes  require  D  to  be  fractional  to  obtain  consistent  results. 
The  fractional  (fractal)  dimension  D  in  these  cases  is  analogous  to  the  Euclidian  dimension  in 
terms  of  its  relation  to  the  metric  properties  of  a  surface. 

In  his  extension  of  the  measuring  concept  outlined  above,  Mandelbrot  starts  with  a  large 
cube  containing  part,  or  all,  of  a  surface  whose  fractal  dimension  is  to  be  determined.  He  defines 
the  side  of  that  cube  as  the  external  scale,  L,  and  then  subdivides  it  into  smaller  cells  of  side 
length  /  counting  the  number  of  smaller  cells  containing  at  least  one  point  on  the  surface.  If  the 
“surface”  is  a  straight  line,  the  number  of  small  cells  through  which  it  passes  is  proportional  to 
(L//)^.  A  flat  plane  will  intercept  a  number  of  small  cells  that  is  approximately  proportional  to 
(IV/)^.  For  a  solid  “surface,”  that  number  of  cells  is  proportional  to  (L//)3.  In  each  example,  the 
exponent  of  (L//)  is  the  Euclidian  dimension  of  the  surface,  confirming  the  analogy  to  the  usual 
measurement  concepts.  However,  the  exponent  assumes  fractional  values  for  rough,  scaling 
objects.  For  example,  a  line  that  has  many  wiggles,  with  wiggles  upon  those  wiggles,  passes 
through  a  number  (N)  of  small  cells  that  is  approximated  by 

N  =  k(L//)D  ,  (2) 

where  k  is  a  proportionality  constant  and  D,  the  fractal  dimension,  is  larger  than  1.  Similarly,  for 
a  rough  plane  with  an  irregular  surface  over  a  wide  range  of  scales,  D  is  larger  than  2. 

Greenside  et  al.  (1982)  attributed  an  algorithm  for  computing  fractal  dimension  based  on 
Mandelbrot’s  (1975)  definition  outlined  above  to  Takens  (1981).  They  found  the  approach  to  be 
impractical  for  surfaces  of  dimension  larger  than  two  because  of  convergence  problems.  Later  in 
this  report,  applications  of  this  approach  to  two-dimensional  scalar  fields  are  presented.  In  these 
applications,  the  scalar  values  were  rescaled  to  fall  within  the  same  numerical  range  as  the  linear 
dimensions  of  the  area  over  which  they  were  measured.  This  allowed  the  use  of  cubes  to  cover 
the  scalar  surface. 

When  the  definition  provided  by  Eq.  (2)  is  reduced  to  two  dimensions,  a  square  and  smaller 
squares  replace  the  cube  and  subcubes;  the  definition  then  can  be  used  to  examine  the  lines 
formed  when  a  plane  intersects  a  fractal  surface.  It  is  generally  accepted  (Mandelbrot,  1983)  that 
if  a  fractal  shape  is  intersected  by  a  flat  surface  of  lower  dimension,  the  intersection  is  also  a 
fractal.  The  dimension  of  that  fractal  differs  from  that  of  the  original  object  by  the  integer  differ¬ 
ence  between  the  dimensions  of  the  two  spaces.  This  fact  could  be  used  to  examine  the  distribu¬ 
tions  of  scalars  in  space  by  applying  it  to  an  isoline  (line  of  constant  scalar  value),  which  is  the 
intersection  of  a  surface  whose  height  is  proportional  to  the  value  of  some  parameter  at  that  x,  y 
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location  in  the  underlying  plane.  Thus,  the  cell-counting  method  could  be  used  to  infer  fractal 
dimension  of  a  parameter  distribution  from  the  fractal  dimension  of  an  isoline.  There  are  pitfalls 
in  this  method  when  the  intersecting  plane  is  not  one  of  constant  parameter  value,  or  the  units  of 
X  and  y  are  different. 

The  cell-counting  concept  is  reasonably  direct  and  easily  applied  to  an  array  of  discrete 
numbers  (e.g.,  a  digitized  image).  The  constants  for  the  logarithmic  form  of  Eq.  (2)  are  found  by 
linear  regression: 

log  (N)  =  log  (k)  +  D  log(  U/)  (3) 

In  practice,  we  want  as  many  values  for  N  as  possible.  For  a  discrete  array,  each  value  must  cor¬ 
respond  to  an  integer  value  of  (  L//).  Therefore,  the  calculations  that  are  described  later  used  a 
large  area  that  was  defined  so  that  L  had  many  divisors.  The  large  number  of  divisors  was 
obtained  hy  defining  L  as  a  product  of  powers  of  2,  3,  and  5. 

2.  Area/Perimeter  Methods 

Another  method  for  calculating  fractal  dimension  that  is  related  to  the  space-fillingness  of  a 
curve  uses  the  relationship  between  the  area  A  enclosed  by  some  shape  and  its  perimeter  P.  A 
snuxnh  perimeter  encloses  a  larger  area  than  the  same  wiggly  length.  For  smooth  curves,  the 
perimeter  P  is  proportional  to  the  square  root  of  the  enclosed  area,  but  a  fractal  shape  encloses  an 
area  related  to  its  perimeter  by 


where  D  is  the  fractal  dimension.  Eqn.  (4)  applies  to  isoplcths  and  their  enclosed  areas,  so  it  can 
also  be  used  for  estimating  the  fractal  dimension  of  a  scalar  field.  The  data  available  for  this 
study  did  not  lend  itself  easily  to  the  identification  of  isopleths  on  a  fine  scale;  thus,  this  tech¬ 
nique  was  not  implemented. 

3,  Scaling  and  Probability  Distributions 

Although  fractals  are  geometric  shapes,  their  potential  applications  as  descriptors  of  natural 
phenomena  frequently  arise  in  connection  with  highly  irregular  spatial  or  temporal  distributions 
of  some  nongeometric  quantity  such  as  kinetic  energy,  aerosol  concentration,  or  temperature. 
Jones  et  al.  ( 1991 )  distinguish  between  two  distinct  types  of  fractal  dimension: 

•  The  fractal  dimension  of  the  graph  of  some  fluctuating  function 

•  The  fractal  dimension  of  the  set  of  points  that  define  active  regions  of  fluctuation. 

The  .second,  geometric  definition  is  addressed  by  the  two  methods  described  in  the  preceding 
sections.  While  it  is  possible  to  convert  a  problem  involving  the  first  type  of  fractal  dimension 
into  one  involving  the  second  type  by  examining  shapes  of  isopleths  as  described  above,  it  is 
often  more  useful  to  invoke  statistical  descriptions.  For  example,  one  common  statistical 
de.scriptor  is  the  probability  distribution  Pr(Ac(Ar)  >  q]  associated  with  differences  in  parameter 
value  over  a  specified  distance.  Pr  (r  >  s  ]  is  the  probability  that  the  argument  inequality  is  true; 


7 


Ac(x)  is  the  absolute  value  of  the  difference  in  c  at  two  points  separated  by  a  distance  x.  These 
probability  distributions  vary  according  to  the  separation  distance.  Generally,  the  likelihood  of 
exceeding  some  specified  difference  increases  with  increasing  separation  between  points.  For 
many  natural  phenomena,  the  probability  distributions  retain  the  same  functional  form  for  differ¬ 
ent  separations,  but  exhibit  a  change  of  scale.  This  is  because  many  natural  fields  vary  in  such  a 
way  that  the  fluctuations  look  qualitatively  the  same,  regardless  of  magnification.  Large-scale 
fluctuations  in  a  scaling  field  are  qualitatively  the  same  as  the  middle-scale  fluctuations  embed¬ 
ded  within  them,  which  are  in  turn  similar  to  the  still  smaller  fluctuations. 

When  a  field  is  scaling,  the  probability  distributions  for  different  separations  are  related  as 
follows: 

Pr  [Ac(XAr)  >  q]  =  Pr  [  >  ^Ac(Ar)  >  q]  .  (5) 

Equation  5  shows  that  the  fluctuations  in  c  over  small  separation  distances  Ar  are  related  to 
those  at  larger  separations  (X  Ar)  by  the  factor  The  exponent  H  is  the  Hausdorff  dimension, 
which  is  related  to  the  fractal  dimension  D  and  the  number  of  Euclidian  dimensions  E  of  the 
space  where  the  field  is  plotted  (Voss,  1988)  by: 

D  =  1  +  E  -  H  .  (6) 

For  example,  a  time  series  (where  time  is  the  dimension  against  which  the  variable  is  plotted 
and  E  =  1)  has  a  fractal  dimension  of  2-H.  When  the  probability  of  the  difference  exceeding 
some  selected  value  is  directly  proportional  to  the  separation,  as  it  is  fora  steadily  increasing  or 
decreasing  variable,  then  H  is  1  and  the  fractal  dimension  is  also  1^ — that  is,  a  smooth  line  for  the 
time  series.  If  H  approaches  0,  the  fractal  dimension  approaches  2  and  the  difference  in  value  at 
two  points  is  almost  as  likely  to  exceed  a  specified  amount  for  closely  spaced  points  as  it  is  for 
widely  spaced  points.  This  is  characteristic  of  a  very  rough,  space-filling  time-series  graph. 

Schertzer  and  Lovejoy  (1983)  assumed  a  functional  form  for  the  upper  tails  of  atmospheric 
probability  distributions  in  order  to  use  Eqs.  (5)  and  (6)  to  estimate  fractal  dimension  from 
observed  data.  Often,  the  tail  of  a  distribution  is  most  interesting  because  it  involves  the  large, 
infrequent  fluctuations.  They  assumed  that  the  upper  tail  had  a  hyperbolic  form;  that  is,  for 
large  q: 

Pr[Ac>q)  =  F(q)  (7) 

where  Schertzer  and  Lovejoy  (1983)  define  F(q)  by: 

F(q)  =  Kq-«  .  (8) 

The  constants  H  and  a  are  estimated  from  the  data  after  substituting  Eq.(8)  in  Eq.  (5)  and  taking 
logarithms: 

log  {  Pr  (Ac  (X.Ar)  >  q]  )  =  k  -i-  H  log(A.)  -  a  log(q)  .  (9) 

This  equation  is  applied  to  observed  data  as  follows: 

n)  a  and  [k  +  H  log(L)l  are  determined  by  linear  regression  from  the  upper  tail  of  the  dis¬ 
tribution  for  each  X. 
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(2)  The  average  a  is  determined;  if  individual  values  of  a  do  not  differ  widely,  the  average 
a  is  used  to  get  an  equation  relating  X  and  H  for  each  value  of  q  so  that  linear  regres¬ 
sion  can  be  used  to  determine  H. 

(3)  D  is  then  obtained  from  Eq.  (6). 

The  hyperbolic  functional  form  for  the  upper  tail  was  chosen  by  Schertzer  and  Lovejoy 
(1983)  because  it  appears  to  fit  atmospheric  distributions  that  they  analyzed.  Other  functional 
forms  can  be  used  in  much  the  same  way  as  described  above.  The  important  point  is  that  the 
probability  distributions  are  the  same  if  the  separation  distances  and  parameter  values  are  scaled 
as  in  Eq.  (5). 

4.  Spectral  Analysis 

The  method  just  described,  which  relates  the  difference  in  the  value  of  a  parameter  at  two 
points  to  the  separation  between  those  points,  is  closely  related  to  the  estimation  of  spatial  auto¬ 
correlation.  Not  surprisingly,  therefore,  spectral-analysis  methods  can  be  used  to  estimate  the 
Hausdorff  and  fractal  dimensions  for  fractal  Brownian  functions  defined  by  Eq.  (5).  For  exam¬ 
ple,  the  exponent  of  the  power  spectrum  (-y )  for  a  scaling  time-series  graph  is  related  to  the 
Hausdorff  dimension  (Saupe,  1988)  by 

y  =  2H-hE  (10) 

When  we  combine  Eqs.  (5)  and  (10)  and  substitute  E  =  1  (for  a  time-series  graph),  we  get 
y  =  5  -  2  D.  The  spectral  density  P(f)  of  a  scaling  variable  is  proportional  to  where  f  is 
frequency. 

The  spectral  relationship  can  also  be  used  for  estimating  the  fractal  dimension  of  a  scaling, 
isotropic  field  in  two  or  three  dimensions.  Fast-Fourier-transform  (FFT)  methods  give  the  ampli¬ 
tudes  of  the  Fourier  terms  in  wave-number  space.  The  spectrum  is  then  fit  without  regard  to 
direction.  In  essence,  the  procedure  averages  the  amplitudes  around  circles  (or  spheres)  of  con¬ 
stant  wave-number  radius.  Pentland  (1984)  used  the  spectral  method  to  characterize  the  fractal 
dimensions  of  16  x  16  pixel  subsections  of  photographs,  and  used  the  results  as  a  basis  for  dis¬ 
tinguishing  between  different  picture  areas  according  to  their  texture.  For  these  two-dimensional 
images,  the  power  density  scaling  is  proportional  to 

The  spectral  approach  proved  the  least  successful  of  those  that  were  implemented  for  these 
studies.  However,  an  inverse  Fourier  transform  approach  described  later  proved  quite  useful  for 
generating  test  arrays  of  known  fractal  dimension. 


5.  Multiresolution  Feature  Analysis 

If  fractals  represent  a  composite  of  similar  “features”  over  a  broad  range  of  scales,  then  it 
should  be  possible  to  use  this  fact  to  estimate  fractal  dimension.  This  is  the  essence  of  the 
approach  described  by  Jones  et  al.  (1991).  The  following  description  of  the  methodology,  taken 
from  Jones  et  al.  (1991),  is  a  variation  of  the  method  based  on  probability  distributions  that  was 
discussed  earlier.  Instead  of  simply  using  changes  in  the  value  of  the  dependent  variable  as  a 
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function  of  distance,  they  express  those  changes  in  the  form  of  a  set  of  differences  that  can  be 
used  to  define  geometric  “features”  of  the  distribution.  They  also  degrade  the  image  resolution 
to  get  different  separation  distances.  This  approach  provides  a  method  by  which  fractal  dimen¬ 
sion  can  be  more  easily  related  to  specific  physical  characteristics  of  the  field. 


The  process  begins  by  producing  an  image  of  half  the  original  resolution  through  low-pass 
filtering.  Correlation  with  a  centered  3x3  matrix  gives  a  value  for  each  point.  These  new  val¬ 
ues  are  resampled  at  every  other  grid  point  in  both  directions,  giving  one  point  for  each  four  in 
the  original  array  and  a  new  image  of  half  the  resolution  (and  about  half  the  grid  points  in  each 
direction).  The  process  is  repeated  to  obtain  a  set  of  images  covering  a  range  of  resolutions. 

Jones  et  al.  (1991)  found  empirically  that  preprocessing  of  the  original  image  to  introduce 
slight  smoothing  helped  them  obtain  better  scaling  in  their  analyses  of  LANDS  AT  imagery.  Fol¬ 
lowing  their  example,  the  data  were  initially  smoothed  with  the  following  matrix: 


L  = 


0 


0.031  0.044  0.031 
0.044  0.7  0.044 

0.031  0.044  0.031 


(11) 


The  smoothings  to  obtain  the  reduced-resolution  imagery  used  the  matrix 


L 


1 


0.0625  0.125  0.0625* 
0.125  0.25  0.125 

0.0625  0.125  0.0625 


(12) 


After  the  set  of  reduced  resolution  images  have  been  obtained,  a  high-pass  filter  is  applied 
to  each  image.  The  high-pass  filter  can  be  interpreted  as  a  “feature  detector”:  It  returns  large 
positive  or  negative  numbers  when  superimposed  on  a  pattern  in  the  image  that  corresponds  to 
the  feature  it  is  designed  to  detect.  For  example,  the  following  matrices  detect  “edges”  oriented 
vertically  and  peaks  (or  holes)  in  the  pattern: 


and 


H 


vert  edge 


10  -1 
10  -1 
10  -1 


1-1-1 
1  8-1 
1  -1  -1 


(13) 


(14) 


The  matrices  can  be  combined.  For  instance,  the  outputs  from  four  edge  detectors  (for  vertical, 
horizontal,  and  two  orientations  of  diagonal  edges)  can  be  summed  to  find  regions  of  strong 
gradients. 

After  arrays  of  feature  intensities  have  been  obtained  for  the  different  resolution  images,  the 
maxima  and  minima  are  identified.  Then  the  numbers  of  these  extremes  whose  absolute  values 
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exceed  various  thresholds  are  counted.  Plots  of  these  thresholds  and  numbers  of  extrema  exceed¬ 
ing  them  for  each  resolution  can  serve  as  a  basis  for  estimating  fractal  dimension.  The  fractal 
dimension  of  the  graph  of  a  fluctuating  function  C  is  the  measure  that  is  most  closely  related  to 
the  scaling  properties  of  C.  According  to  Eq.  (5)  self-similar  scaling  processes  display  amplitude 
fluctuations  (AC)  whose  magnitudes  tend  to  be  proportional  to  the  ratio  of  their  spatial  scales  \ 
raised  to  a  power,  H.  Jones  et  al.  (1991)  refer  to  H  as  the  similarity  parameter.  In  essence,  the 
process  described  above  identifies  the  number  (or  probability)  of  fluctuations  of  certain  type 
exceeding  various  thresholds  at  different  resolutions. 

The  process  described  above  provides  an  estimate  of  (nq)i ;  the  number  of  exceedances  (for 
a  given  total  area)  of  threshold  q  for  the  i*  resolution  XiArp,  where  atq  is  the  cell  size  in  the  origi¬ 
nal  data  set.  If  the  distribution  is  self  similar,  then  the  number  of  features  exceeding  a  given 
threshold  q  will  satisfy  a  relationship  of  the  following  form: 

^2  -H 

(nq).  =  q\.  ,  (15) 

Jones  et  al.  (1991)  choose  a  series  of  values  for  H  and  estimate  fractal  dimension  based  on 
which  of  those  values  produces  the  best  coincidence  of  the  graphs  for  the  different  resolutions. 
Examples  of  the  method  are  given  later. 

6.  Anisotropic  Effects 

The  above  discussions  all  assume  that  the  relationships  are  both  independent  of  direction 
and  scale.  This  is  not  necessarily  the  case.  Differences  between  scaling  of  feature  dimensions  in 
the  vertical  and  horizontal  directions  seem  likely  to  occur  in  a  medium  that  is  as  highly  stratified 
as  the  atmosphere.  It  turns  out  that  many  of  the  fractal  properties  of  the  atmosphere  arise  from 
the  cascade  of  energy  down  through  the  various  scales  of  motion.  A  fractal  dimension  that 
applies  over  the  complete  range  of  scales  implies  that  the  ratios  of  the  number  of  eddies  to 
subeddies  remains  constant  over  that  range,  which  need  not  be  true.  While  this  report  has  not 
attempted  to  incorporate  these  effects  into  the  analyses  presented  later,  completeness  requites 
that  they  be  discussed. 

a.  Self-Affine  Distributions 

The  analytical  techniques  discussed  above  are  appropriate  for  shapes  with  self-similar  char¬ 
acteristics,  that  is,  for  those  cases  where  the  irregularities  at  all  scales  are  similar  in  the  Euclidian 
geometry  sense.  The  smaller-scale  features  are  contracted  (and  perhaps  rotated)  versions  of  the 
larger  irregularities,  but  they  maintain  their  proportions  and  the  corresponding  angles  do  not 
change.  These  cases  are  a  subset  of  a  more  general  class  of  shapes  that  occur  naturally,  the  self- 
affine  fractals,  where  the  smaller-scale  irregularities  are  affine  transformations  of  those  at  larger 
scales.  Barnsley  and  Sloan  (1988)  described  affine  transformations  as  “combinations  of  rotations, 
scalings,  and  translations  of  the  coordinate  axes  in  n-dimensional  space.”  Voss  (1988)  defines 
self-affinity  as  nonuniform  scaling,  i.e.,  different  scaling  along  different  axes  (directions). 
Mandelbrot  (1985)  pointed  out  that  the  methods  described  above  for  calculating  fractal  dimen¬ 
sion  are  not  always  appropriate  for  this  more  general  class.  In  particular,  he  notes  that  the  time 
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history  of  a  scalar  involves  a  choice  of  units  for  both  time  and  the  scalar,  different  results  may  be 
obtained  for  different  choices  of  units.  There  will  be  local  and  global  fractal  dimensions.  In  fact, 
Mandelbrot  (1985)  defines  three  versions  of  fractal  dimension  corresponding  to  different  meth¬ 
ods  for  estimating  them.  The  three  different  dimensions  are  the  same  for  self-similar  sets,  but 
differ  for  self-affine  sets.  Mandelbrot  (1985)  explained  that  (for  the  time  series  of  some  param¬ 
eter,  B) 

../'square,"  "distance,"  and  "circle”  are  vital  notions  in  isotropic  geometry, 
but  they  are  meaningless  in  affine  geometry.  More  precisely,  they  are  meaning- 
fill  for  relief  cross-sections,  but  are  meaningless  for  [B]  noises,  because  the 
units  along  the  t-axis  and  the  B-axis  are  set  up  independently  of  each  other, 
hence  At  and  aB  cannot  be  combined.  There  is  no  intrinsic  meaning  to  the 
notion  of  equal  height  and  width,  a  square  cannot  be  defined.  Similarly,  a  cir¬ 
cle  cannot  be  defined  because  its  square  radius  -  At^  +  aB^  would  have  to 
combine  the  units  along  both  axes.  Furthermore,  one  cannot ’walk  a  compass’ 
along  a  self-affine  curve  because  the  distance  covered  by  each  step  combines  a 
At  and  a  aB. 

The  above  discussion  emphasizes  the  care  that  must  be  taken  when  applying  what  are  essen¬ 
tially  geometric  measurement  techniques  to  the  problem  of  estimating  a  fractal  dimension. 
Although  two  of  the  techniques  described  by  Mandelbrot  (1985)  give  results  in  the  small  scale 
that  are  consistent  with  the  results  obtained  for  self-similar  sets,  it  is  wise  to  apply  the  techniques 
carefully  to  avoid  the  inconsistencies  enumerated  above.  Even  then,  anisotropy  in  the  processes 
producing  the  sets  can  make  new  measures  necessary.  One  of  these  is  discussed  in  the  next 
subsection. 

b.  Elliptical  Dimension 

The  horizontal  scales  of  atmospheric  motion  are  much  less  constricted  than  the  vertical. 

This  causes  anisotropy  in  atmospheric  turbulence,  especially  at  the  larger  scales.  It  also  suggests 
that  the  transfer  of  energy  from  large  eddies  to  smaller  eddies  is  not  likely  to  be  self-similar,  but 
rather  self-affine,  even  though  the  vertical  and  horizontal  distance  units  are  the  same.  Schertzer 
and  Lovejoy  (1983,  1987;  also  Lovejoy  and  Schertzer,  1986)  noted  the  following  important  facts 
regarding  mesoscale  atmospheric  processes: 

•  The  energy  spectrum  is  scaling  (i.e.,  it  has  the  form  k^h,  where  k  is  the  wave  number 
and  ph  =  5/3  is  the  appropriate  value  for  wave  numbers  in  the  horizontal  plane). 

•  The  energy  spectrum  for  wave  numbers  in  the  vertical  plane  is  also  scaling,  but 
anisotropy  makes  the  relevant  exponent,  pv»  quite  different:  pv  =  1 1/5. 

•  There  is  extreme  variability,  because  aaive  regions  that  account  for  most  of  the  energy 
and  moisture  flux  are  sparsely  distributed. 

To  deal  with  the  above  observations,  Lovejoy  and  Schertzer  (1986)  modified  the  self-simi¬ 
lar  concepts  to  apply  them  to  the  anisotropic  case.  In  extending  the  concepts,  they  assumed  that 
there  is  a  constant  energy  flux  over  the  range  of  scales  of  interest,  and  that  there  are  rules  for 
describing  how  the  statistical  properties  of  eddies  are  transformed  from  one  scale  to  another.  In 
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dealing  with  anisotropic  eddies,  both  shape  and  dimension  must  be  characterized.  Schertzer  and 
Lovejoy  (1987)  represented  the  anisotropic  scaling  in  a  manner  similar  to  that  discussed  earlier 
for  the  isotropic  case,  except  that  horizontal  and  vertical  scaling  were  considered  separately,  that 
is. 


Pr  ( Ac(XAx)  >  q  ]  =  Pr  [  Ac(Ax)  >  q  ] 

(16a) 

Pr  [  Ac(XAy)  >  q  ]  =  Pr  [  Ac(Ay)  >  q  ] 

(16b) 

Pr  1  Ac(XAz)  >  q  ]  =  Pr  [  Ac(Az)  >  q  ] 

(16c) 

Pr  [  Ac(T  (Ar)  >  q  ]  =  Pr  [  X”*'  Ac(Ar)  >  q  ] 

(17) 

where 

X  0  0 
0X0 

0  0  x”* 

Hz  =  (Hh/Hv)  (19) 

and 


The  matrix  T  produces  a  magnification  overall,  with  stretching  in  the  z  direction;  it  trans¬ 
forms  the  probability  distributions  and  introduces  an  elliptical  geometry  to  account  for  the  differ¬ 
ent  horizontal  and  vertical  scalings.  Figure  1  shows  how  the  magnification  and  stretching  relates 
the  small,  vertically  oriented  shapes  (eddies  in  this  case)  to  a  large,  horizontally  oriented  shape. 
The  transformation  changes  the  volume  by  a  factor  ,  where  Del  =  2  +  Hz.  Lovejoy 

and  Schertzer  (1986)  called  Del  the  elliptical  dimension.  In  an  isotropic  system,  =  1 ,  and 
therefore  the  elliptical  dimension  equals  3.  For  two-dimensionally  isotropic  sets  Hz  =  0  and  Dd 
=  2.  Schertzer  and  Lovejoy  (1987)  described  a  technique  called  functional  box  counting  for 
determining  elliptical  fractal  dimension  from  observed  data. 

C.  HOW  FRACTALS  MAY  DEVELOP 

To  this  point,  the  discussion  has  focused  on  definitions  of  fractals  and  related  concepts, 
while  giving  little  attention  to  their  underlying  physical  causes.  At  least  three  types  of  evidence 
indicate  a  connection  between  fractals  and  atmospheric  processes: 

•  The  structure  of  turbulent  eddies  observed  in  flow  visualization  experiments 
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Similarity  arguments  from  classical  turbulence  theory 
Detailed  numerical  flow  simulations. 


These  are  discussed  in  the  following  sections. 

1.  Observational  Evidence 

a.  Shear<Layer  Flow 

Van  Dyke  (1982)  showed  several  examples  of  structure  in  the  turbulent  layer  that  forms  in 
the  shear  zone  between  fluids  flowing  at  different  velocities.  The  structure  that  commonly  occurs 
in  such  flows  is  evident.  Reynolds  (1985)  provided  the  basis  for  a  qualitative  discussion  of  the 
processes  involved.  Although  his  notes  do  not  mention  fractals,  they  do  give  evidence  of  scaling 
and  the  development  of  self-affine  structures  of  different  scales  in  the  motion  field. 

Figure  2  is  based  on  Reynolds’  (1985)  schematic  diagram  showing  free  shear  flow.  Initially, 
vorticity  is  produced  at  the  tip  of  the  separator  that  is  used  to  generate  the  flow.  An  unstable  two- 
dimensional  shear  layer  is  formed.  Instabilities  excited  by  slight  vibrations  grow  rapidly  so  that 
vorticity  is  quickly  concentrated  in  nearly  discrete  vortices,  oriented  generally  perpendicular  to 
the  flow  direction.  Not  all  vorticity  is  in  the  major  features;  some  remains  in  “braids”  connecting 
the  major  vortices  (Figure  2),  where  it  is  stretched  as  the  large-scale  vortices  “wind  in”  the  fluid 
between  them.  This  intensifies  the  vorticity  in  the  braids,  forming  new  vortices  with  axes  aligned 
along  the  principal  strain  direction,  as  shown  in  Figure  3.  These  vortices  undergo  the  same  pro¬ 
cesses  described  above,  but  on  a  smaller  scale  and  with  different  orientation.  One  can  expect  that 
the  new,  smaller  vortices  also  contain  irregularities  where  new  minibraids  form  and  stretch  vor¬ 
tices  along  new  strain  axes.  These  new  structures  can  undergo  similar  deformation  on  a  still 
smaller  scale,  so  that  the  processes  continue  to  produce  self-affine  structures  from  the  outer  scale 
defined  by  the  largest  vortices  down  to  the  scale  of  viscous  dissipation. 

The  consequence  of  this  process  appears  to  be  a  set  of  vortices  on  vortices  on  vortices  and 
so  on,  much  like  Richardson’s  (1922)  description.  Analogous  reasoning  can  be  applied  to  other 
flow  types  such  as  jets,  boundary  layers,  and  wakes,  which  are  also  characterized  by  regions  of 
strong  shear  and  concentrated  vorticity. 

b.  Turbulent  Diffusion 

Procaccia  (1984)  reviewed  the  effects  of  fractal  turbulence  structures  on  turbulent  diffusion, 
fluctuations  of  passive  scalars,  electromagnetic  wave  propagation,  and  cloud  perimeters.  He 
examined  the  behavior  of  the  interparticle  separation  distance,  R  =  -  ri,  of  two  points  caused 

by  their  relative  velocity,  V(t)  =  Vj  -  V2‘. 

t 

R(t)  =  R(0)  +  J  VfTldt  .  (21) 

0 
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Entrainment  of 
external  fluid 


FIGURE  2  VORTEX  BRAIDS  BETWEEN  SHEAR  FLOW  EDDIES 
After  Reynolds,  1985 
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FIGURE  3  ENERGY  TRANSFER  TO  SMALLER  EDDIES  BY  VORTEX 
STRETCHING 

After  Reynolds.  1985 
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For  isotropic  turbu.c/ice,  the  ensemble  average  separation — represented  by  angle  brackets,  <>  is 
constant  and  equal  to  the  initial  separation  because  <  V  (t )  >  =  0 ,  but  variance,  <  R2>, 
changes: 


(22) 


Procaccia  (1984)  asserted  that,  although  the  correlation  <  V  (t )  •  V  (x  )>  is  nonstationary, 
some  function,  g(x),  of  scaled  time  variables  exists  such  that 


<  V(t)-V(t)>=(  V(t)-V(t)>g 


(23) 


where  tR  is  the  typical  decay  time  for  velocity  differences  over  a  length  scale,  R.  Substituting  in 
Eq.  (22)  provides  asymptotic  predictions  at  extreme  times: 


s  <'■>-{ 


<  V(t).V(t))  t 

<  V(t).V(t)) 


(24) 


The  diffusivity,  d  <  R2>  /dt,  can  be  determined  from  <  V(t)*V(t)>  when  R  is  in  the  inertial 
range.  The  ratio  of  the  separation  distance,  R,  to  a  typical  velocity  difference  over  that  distance 
defines  tR.  For  the  “homogeneous  fractal  model”  of  turbulence  described  below. 


<v.v)  -(eR)^  (7  )  ^  • 


(25) 


where  /q  is  the  outer  scale  of  the  turbulence.  Procaccia  (1984)  used  this  r,ssumption  to  obtain 


jd 

dt 


(26) 


where  R  is  the  root-mean-square  separation.  For  space-filling  turbulence,  where  D  =  3,  the 
above  expressions  reduce  to  the  classical  “4/3  law.” 

Hentschel  and  Procaccia  (1983)  examined  Gifford’s  (1957)  and  Richardson’s  (1926)  two- 
point  diffusion  data  and  obtained  estimates  for  D  between  2.5  and  2.78,  from  which  they  con¬ 
cluded  that  these  were  reasonable  values  that  supported  the  “fractaily  homogeneous  turbulence” 
concept. 
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c.  Cloud  and  Rain  Perimeters 

As  noted  earlier,  Lovejoy  (1982)  found  that  cloud  and  rain  area  perimeters  were  fractals 
over  a  range  of  six  orders  of  magnitude — from  about  1  km^  to  10^  km^.  The  fractal  dimension 
of  the  perimeter  was  D  =  1.35  ±  0.05,  so  the  fractal  dimension  of  the  cloud  surface  was  2.35  ± 
0.05  in  the  isotropic  case.  Procaccia  (1984)  considered  how  a  surface  defining  the  outer  boun¬ 
dary  of  the  cloud  would  be  distorted  by  turbulent  diffusion,  and  concluded  that  the  fractal 
dimension,  Dc,  of  the  cloud  perimeter  is  given  by 

Dc  =  (ll-D)/6  ,  (27) 

where  D  is  the  fractal  dimension  of  the  turbulence.  Usin_g  a  value  of  2.6  for  the  fractal  dimension 
of  the  turbulence  (Hentschel  and  Procaccia,  1983)  gives  D  =  1.4,  in  agreement  with  Lovejoy’s 
(1982)  value. 

d.  The  Boundary  of  a  Diffusing  Scalar  In  a  Turbulent  Flow 

Prasad  and  Sreenivasan  (1990)  experimentally  examined  several  aspects  of  turbulence  for 
fractal  characteristics.  They  measured  the  concentration  field  of  a  diffusing  fluorescent  material 
in  water.  Their  technique  used  a  sheet  of  laser  light  to  induce  fluorescence.  The  fluid  was 
scanned  very  rapidly  and  two-dimensional  images  obtained  from  closely  spaced  planes.  Subse¬ 
quent  computer  processing  provided  the  three-dimensional  distribu'ion  of  the  material.  Among 
the  characteristics  that  Prasad  and  Sreenivasan  (1990)  examined  was  the  shape  of  the  interface 
between  the  fluorescent  material  and  the  surrounding  water.  They  found  that  the  interface  was 
very  rough,  and  generally  self-similar  with  a  fractal  dimension  of  2.35  ±  0.04. 

e.  Elliptical  Fractals 

Schertzer  and  Lovejoy  (1983)  introduced  the  concept  of  elliptical  fractal  dimension  because 
it  provides  a  smooth  transition  from  the  very  large,  horizontally  oriented,  “Hadley-like”  cells 
down  to  the  vertically  oriented  convective  cells.  An  analysis  (Lovejoy  et  al.,  1987)  of  radar  pre¬ 
cipitation  data  showed  that  rainfall  is  distributed  in  space  with  an  elliptical  fractal  dimension  of 
about  2.22.  The  elliptical  fractal  dimension  represents  the  ratio  of  horizontal  contraction  to  verti¬ 
cal.  so  this  suggests  that  rainfall  is  more  horizontally  stratified  than  are  the  air  motions,  which  the 
authors  estimated  as  having  a  fractal  dimension  of  about  2.56  (Lovejoy  et  al.,  1987). 

2.  Classical  Turbulence  Theory 

The  description  by  Frisch  et  al.  (1978)  of  the  “p-theory”  of  turbulence  provides  a  good 
physical  picture  of  the  turbulent  cascade.  The  energy  spectrum  E(k)  is  defined  as  the  kinetic 
energy  per  unit  mass  per  unit  wave  number  k.  For  purposes  of  argument,  consider  a  spectrum  of 
eddies  beginning  at  the  largest  scale,  /o>  where  the  energy  is  introduced,  and  proceeding  to  suc¬ 
cessively  smaller,  discrete  sizes,  /n,  as  follows: 

kn=  l//„  =  2"//o  n=0,  1,2,3,...  (28) 
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If  En  is  the  discretized  kinetic  energy  per  unit  mass  in  the  wave  numbers  near  kn,  it  can  be 
defined  in  terms  of  a  characteristic  velocity  Vr.  This  characteristic  velocity  is  not  total  velocity, 
but  the  velocity  difference  for  the  eddy  size  /n- 

E_  ~v^ 

^  "  .  (29) 

The  eddy  turnover  time  tn  is  given  by 

in~^n/Vn  ,  (30) 

Frisch  et  al.  (1978)  defined  an  energy  (per  unit  mass)  transferrate,  En,  from  eddies  of  wave 
number  kn  to  kn+i  for  the  inertial  subrange,  where  tn  is  much  greater  than  the  viscous  dissipation 
time,  /n^/v  (V  is  kinematic  viscosity),  and  much  less  than  the  characteristic  time  for  the  larger 
scale  motions, 

En  —  En/tn  ”*  (Vn)V^ n  •  (31) 

For  statistically  stationary  turbulence,  energy  introduced  at  large  scales  (/q)  transfers  to  suc¬ 
cessively  smaller  scales  until  the  dissipation  scale,  /j,  is  reached  and  the  energy  dissipation  rate, 
En,  is  a  constant,  t.  Solving  for  Vn  and  En  gives 


(e/„) 


2 

E„- 


Equations  (33)  and  (34)  are  the  same  as  Kolmogorov’s  results.  A  Fourier  transformation 
provides  the  wave-number  spectrum: 

1  _  1 

E(k)~(E)^k  ^  .  (34) 

Equation  (33)  can  be  rewritten  to  show  that  the  energy  per  unit  mass  per  unit  volume  scales 
according  to  (///o  by  noting  that 

=  (35) 

Then, 

111  1 

=  ,  (36) 

There  is  a  tacit  (but  very  important)  assumption  in  the  above  relationship:  All  eddy  sizes 
are  assumed  to  be  spread  more  or  less  uniformly  throughout  the  same  volume.  For  greater 
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generality,  Frisch  et  al.  (1978)  assumed  that  the  smaller  eddies  are  less  space-filling  than  the 
larger  ones.  The  behavior  of  free  shear  layer  turbulence  described  in  the  preceding  section  quali¬ 
tatively  supports  this  assumption.  Frisch  et  al.  (1978)  define  a  parameter  to  characterize  the 
degree  of  space-fillingness.  For  a  wholly  space-filling  cascade  of  the  type  discussed  above,  the 
number  of  eddies  of  size  /n  =  is  2^  times  as  those  of  size  /n+i  because  eight  eddies  of 
dimension  //2  fill  the  volume  occupied  by  one  of  dimension  /  Frisch  et  al.  (1978)  defined  p  to 
be  the  average  fraction  of  the  /n  volume  filled  by  /n+i  eddies.  Eddies  of  size  /n  fill  only  a  frac¬ 
tion,  Pn.  of  the  total  space  occupied  by  eddies  of  size  /q. 

Pn  =  (N/23)n  ,  (37) 

where  N  is  the  average  number  of  eddies  formed  by  each  eddy  of  the  preceding  generation,  and 
of  course  N  <  8. 

Frisch  et  al.  (1978)  assumed  that  eddies  of  the  (n+l)!^  generation  are  positionally  correlated 
with  those  of  the  n^  generation  by  embedding  or  attachment,  so  that  the  region  near  where  an 
eddy  is  formed  becomes  an  “active  region”  for  the  cascade  to  smaller  sizes.  The  evidence  from 
mixing-layer  turbulence  suggests  that  attachment  is  more  likely  than  embedding.  Frisch  et  al. 
(1978)  also  tacitly  assumed  that  the  average  number  of  eddies  formed  (N)  does  not  vary  system¬ 
atically  according  to  eddy  size.  This  is  a  weakness  of  this  “3  model”  of  turbulence. 

If  the  kinetic  energy  per  unit  mass  associated  with  scales  on  the  order  of  /n  is  redefined  in 
terms  of  active  regions  only,  then 


(38) 


The  characteristic  energy  transfer  time  is  still  /n/vn  if  the  smaller  eddies  arise  from  the 
internal  dynamics  of  the  larger  eddies  that  produce  them,  but  the  steady-state  assumption  leads  to 
an  adjusted  rate  of  energy  transfer  in  the  inertial  range;  that  is. 


E  B  v^ 

r  n  n  -= 

t - £ 


t_ 


(39) 


The  following  relationships  are  obtained  from  Eq.  (39)  when  we  introduce  D,  defined  by 
N  =2^: 


n  '  n 


(40) 

2/3  (3-D)/3 

(41) 
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and 


-  2/3 -5/3  -(3-D)/3 

E(k)~(er  k  '  (k/^) 

Rewriting  Eq.  (42)  shows  the  scaling  nature  of  the  relationships: 

2/3  (3  -  D)/3  -  2/3  ^33"^  ^33^ 

(/„/V 


(42) 


(43) 


This  is  the  same  as  Eq.  (37),  but  with  a  correction  of  D/3  in  the  exponent  to  account  for  the 
incompleteness  of  space-filling  (intermittency)  in  the  energy  transfer  to  smaller  scales.  The  rela¬ 
tionship  among  N,  p,  and  D  is 

(3  =  (N/23)n  =  (2E>-^)"  .  (44) 

The  above  derivation  demonstrates  the  scaling  properties  of  important  turbulent  parameters. 
It  relates  to  classical  theory  and  observed  physical  factors  that  cause  intermittency.  The  fractal 
dimension  was  introduced  as  a  measure  of  the  assumed  space-filling  properties  of  the  energy 
transfer  process.  When  the  process  is  assumed  to  be  space-filling,  the  fractal  dimension  is  3, 
giving  the  classical  results  of  Kolmogorov. 

Fujisaka  and  Mori  (1979)  estimated  values  for  D  under  the  assumption  that  informational 
entropy  would  be  maximized.  Their  analysis  led  to  an  estimate  of  D  =  2.659  which  is  in  good 
agreement  with  estimates  based  on  observations.  It  also  suggests  that  on  average,  there  are  6.32 
eddies  of  size  /n+i  for  each  eddy  of  size  /n  in  the  above  analysis. 


3.  Numerical  Simulation  of  Small-Scale  Flow  Structures 

Most  fluid  flow  simulations  parameterize  turbulent  effects  and  so  provide  little  evidence  of 
fractal  structures.  However,  Chorin  (1982),  using  vortex  methods  like  those  described  by 
Leonard  (1985),  performed  an  interesting  numerical  experiment  that  bears  a  close  relationship  to 
the  formation  of  vonices  in  the  braids  between  eddies.  He  considered  a  straight  vonex  with  a 
single  perturbation,  much  like  one  of  the  perturbations  on  the  vortex  filaments  in  braids  between 
vortices.  After  a  short  while,  the  vortex  tube  segments  stretched  rapidly.  The  general  orientation 
of  the  pattern  remained  the  same,  but  it  became  contoned  and  complex.  Chorin  (1982)  evaluated 
the  Hausdorff  dimension  of  the  resulting  structure  and  found  it  to  be  on  the  order  of  2.5,  which  is 
consistent  with  values  suggested  by  Mandelbrot  (1977)  for  the  fractal  dimension  of  turbulent 
structures. 
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Ill  IMPORTANT  RESULTS 


A.  SUMMARY  OF  PROJECT  ACCOMPLISHMENTS 

The  research  conducted  during  this  project  has  provided  the  following. 

•  A  coherent  review  and  summary  of  fractal  concepts  and  their  relevance  to  atmospheric 
behavior,  especially  as  they  apply  to  the  turbulent  dispersal  of  aerosols  and  gasses. 

•  Identification  of  several  analysis  techniques  based  on  fractal  concepts  that  can  be 
applied  to  atmospheric  data. 

•  Computer  programs  for  calculating  parameters  necessary  to  estimate  fractal  dimension 
by  methods  described  in  the  literature. 

•  Extension  of  the  multiresolution  feature  analysis  concepts  to  three-dimensional  scalar 
and  vector  fields. 

•  Development  of  computer  codes  for  applying  the  extended  analysis  techniques. 

•  Identification  and  acquisition  of  atmospheric  data  for  analysis. 

•  Analysis  of  some  of  the  available  data. 

•  Development  of  a  methodology  (including  computer  programs)  for  identifying  natural 
patterns  of  variability  in  scalar  and  vector  fields,  so  that  those  patterns  can  serve  as  the 
basis  for  multiresolution  feature  analysis. 

•  Identification  of  promising  approaches  for  future  research. 

The  remainder  of  this  section,  and  the  concluding  sections  provide  examples  of  these 
accomplishments. 

B.  IDENTIFICATION  OF  DATA  SUITABLE  FOR  ANALYSIS 
1.  Background 

The  most  important  purpose  of  this  study  has  been  to  examine  the  spatial  inhomogeneities 
in  plumes  released  in  the  atmosphere.  Obviously,  detailed  observations  of  such  plumes  will 
provide  the  best  possible  data  for  analysis.  It  was  possible  to  obtain  two  such  types  of  data:  One 
data  set  was  collected  with  a  airborne  laser  radar  (lidar)  to  provide  vertical  planar  cross  sections 
of  aerosol  backscatter  from  an  elevated  power  plant  plume,  and  the  other  data  set  (only  identified 
and  obtained  near  the  end  of  the  study)  consists  of  two-dimensional  arrays  of  infrared  transmit¬ 
tance  measurements  through  a  smoke  plume  released  at  the  surface. 

Another  set  of  similar  “data”  was  also  analyzed.  In  this  case,  the  data  were  calculated  to 
have  known  fractal  properties.  While  such  data  will  never  substitute  for  actual  observations,  they 
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have  the  imponant  advantage  of  providing  an  image  that  can  be  used  to  test  the  performance  of 
the  analytical  techniques  that  were  applied  to  the  observations.  The  generation  of  the  test  data 
also  gave  some  insights  into  how  one  might  go  about  generating  artificial  distributions  for  simu¬ 
lating  the  appearance  and  the  effects  of  inhomogeneous  smoke  plumes  on  Army  operations. 

Another  type  of  data  was  obtained  for  analysis,  although  they  were  not  directly  related  to 
inhomogeneities  in  smoke  plumes.  Data  obtained  by  dual  Doppler  radar  observations  of  the 
motions  of  chaff  in  a  convective  atmospheric  boundary  layer  provide  a  detailed  picture  of  those 
atmospheric  motions  that  are  the  ultimate  source  of  the  inhomogeneities  in  the  distributions  of 
interest.  Analyses  of  these  data  have  not  proceeded  as  far  as  those  of  the  scalars,  in  part  because 
they  arc  not  as  directly  relevant  to  project  objectives  as  arc  the  smoke-plume  observations. 

All  the  data  sources  are  discussed  briefly  below.  That  discussion  is  followed  by  a  summary 
of  the  analyses  that  were  completed. 

2.  Scalar  Data 

a.  Random  Brownian  Fractal  Test  Data 

Jones  et  al.  (1991)  and  Saupe  (1988)  describe  a  method  for  generating  spatial  distributions 
of  known  fractal  dimension  by  using  the  spectral  relationships  discussed  earlier  [Eqs.  (5)  and 
(10)].  The  method  is  well-described  in  both  the  sources  cited  above;  Saupe  (1988)  outlines  the 
algorithm  in  a  section  of  computer  pseudocode.  The  description  that  follows  is  based  on  that 
given  by  Jones  et  al.  (1991). 

The  process  begins  by  generating  a  random  array  of  complex  numbers.  The  real  and 
imaginary  parts  each  have  zero  mean  and  unit  variance.  The  examples  shown  in  Figure  4,  are 
based  on  a  256  x  256  array.  The  complex  numbers  in  the  Fourier  plane  are  then  multiplied  by 
their  distance  from  the  origin  (in  wave  numbers)  raised  to  a  power  that  ultimately  determines  the 
fractal  dimension  of  the  distribution  to  be  produced.  If  fi  and  f2  are  the  coordinates  of  a  point  in 
the  Fourier  plane  and  we  wish  to  produce  a  real  number  array  whose  two-dimensional  spectral 
density  S  scales  with  an  exponent  between  0  and  1 ,  i.e., 

S(fi,f2)~(fl2  +  f22)-T^  ,  (45) 

then  the  complex  numbers  in  the  random  array  multiplied  by  the  following  factor,  F  (fl,f2)‘. 

F(fi,f2)  =  (fl2  +  f22)-^'^  »  (46) 

where  from  (Eqs.  (6)  and  (10)J  we  get  the  following  expression  for  the  resulting  fractal 
dimension; 

D  =  (3E  +  2-y)/2  ,  (47) 

or  for  E  =  2,  we  can  choose  y  according  to  the  following  expression  to  obtain  the  desired  fractal 
dimension: 

y=8-2D  .  (48) 
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After  the  complex  numbers  in  the  Fourier  plane  have  been  multiplied  by  the  appropriate 
value  of  F,  the  inverse  Fourier  transform  is  applied  to  obtain  a  new  array  of  complex  numbers  in 
linear  space.  The  real  parts  of  this  array  form  an  array  with  the  desired  fractal  dimension,  at  least 
within  the  limits  of  the  resolution.  Obviously,  there  will  be  a  high-wave-number  cutoff  imposed 
by  the  discrete  nature  of  the  grid. 

Figure  4  shows  three  examples  generated  as  described  above.  Standard  FORTRAN  Sub¬ 
routines  from  Press  et  al.  (1986)  were  used  to  generate  the  desired  random  numbers  and  perform 
the  inverse  fast  Fourier  transform.  The  figure  displays  the  results  as  gray  scale  imagery  with  256 
gray  levels.  The  fractal  dimensions  for  the  three  panels  of  Figure  4  are,  2.3, 2.5,  and  2.7.  The 
same  arrays  used  to  generate  Figure  4  are  analyzed  later  by  three  different  methods  for  calculat¬ 
ing  fractal  dimensions. 

b.  LIdar  Observations  of  Smoke  Plumes 

Uthe  (1983)  described  the  operations  of  the  airborne  lidar  downwind  of  an  Electric  power 
plant  plume  in  Kincaid,  Dlinois.  Figure  5  is  a  schematic  diagram  of  the  operations.  The  aircraft 
flies  at  about  3  km  above  ground  level,  and  the  1 .06-tim  laser  is  pulsed  at  a  rate  that  corresponds 
to  a  distance  of  about  10  m  between  each  vertical  profile  of  aerosol  backscatter.  The  pulse  length 
and  recorder  operation  provide  a  vertical  resolution  of  about  3  m.  The  result  is  an  array  of  values 
corresponding  to  aerosol  backscatter  with  a  spatial  resolution  of  about  3  m  (vertical)  by  10  m 
(horizontal). 

Typically,  the  data  were  collected  in  a  vertical  plane  approximately  normal  to  the  plume  and 
10  to  15  km  downwind  of  the  source.  The  logarithm  of  backscatter  was  recorded  with  8-bit  (256 
units)  resolution.  Before  analysis,  the  data  were  converted  to  a  relative  linear  scale.  The  origi¬ 
nally  recorded  range  of  0  to  255  corresponded  to  16  dB.  Hence,  each  unit  corresponds  to  an 
increase  of  almost  1.5%.  These  data  were  converted  to  a  linear  scale  before  analysis.  Range  cor¬ 
rections  did  not  account  for  the  attenuation  introduced  by  the  relatively  low  aerosol  concentra¬ 
tions  in  the  plume.  This,  and  the  fact  that  aerosol  backscatter  is  not  directly  proportional  to  con¬ 
centration  limit  the  uses  of  the  lidar,  except  for  the  determination  of  geometric  measures.  The 
fractal  analysis  techniques  deal  with  geometric  features  which  should  be  relatively  unaffected  by 
the  lidar’s  limitations  in  estimating  aerosol  concentration. 

The  rectangular  nature  of  the  data  elements,  with  horizontal  dimensions  3.3  x  the  vertical , 
may  have  an  effect  on  the  results.  However,  the  atmosphere  tends  to  horizontally  stratified  with 
more  damped  vertical  motions  and  stronger  vertical  gradients,  so  the  horizontally  stretched  shape 
of  the  data  elements  may  be  more  appropriate  than  a  square  shape. 

Figure  6  shows  three  cross  sections  through  a  smoke  plume  measured  about  1 1  km  down¬ 
wind  of  the  source  between  about  0800  and  0830  on  20  July  1980.  Brighter  regions  indicate 
higher  backscatter.  An  area  of  144  x  144  elements  has  been  selected  from  the  original,  larger 
data  arrays.  Two  of  the  analytical  approaches — spectral  analysis  and  multiresolution  feature 
analysis — are  best  applied  with  data  arrays  whose  dimensions  are  a  power  of  two.  Therefore  128 
X  128  arrays  centered  on  the  plume  were  extracted  from  the  images  shown  in  Figure  6.  The  box¬ 
counting  approach  benefits  by  having  arrays  whose  dimensions  are  divisible  by  many  numbers, 
so  the  full  144  x  144  cells  were  used  for  that  analysis. 
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FIGURE  4  RANDOM  BROWNIAN  FRACTALS 
(a)  FRACTAL  DIMENSION  =  2.3 
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FIGURE  4  RANDOM  BROWNIAN  FRACTALS  (continued) 
(b)  FRACTAL  DIMENSION  =  2.5 


FIGURE  4  RANDOM  BROWNIAN  FRACTALS  (concluded) 
(c)  FRACTAL  DIMENSION  =  2.7 
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FIGURE  5  SCHEMATIC  DIAGRAM  OF  LIDAR  MEASUREMENT  OF  AEROSOL  BACKSCATTER 
FROM  A  SMOKE  PLUME 

c.  Transmissometer  Images 

The  Atmospheric  Transmission  Large-area  Analysis  System  (ATLAS),  described  by 
Bleiweiss  et  al.  (1991),  uses  a  video  imaging  system  to  generate  two-dimensional  arrays  of 
transmittance  values.  The  radiance  of  whatever  is  beyond  the  smoke  plume  being  measured 
serves  as  the  source  for  measurements.  This  means  that  several  conditions  must  be  met  for  the 
data  to  be  reliable.  First,  the  radiance  of  the  background  scene  must  remain  constant  during  the 
measurement  period.  Second,  there  must  be  good  contrast  between  the  background  radiance  and 
that  from  a  wholly  opaque  smoke  plume.  Unlike  the  lidar,  backscatter  data  presented  earlier 
(which  depend  on  concentration),  the  ATLAS  estimates  transmittance  from  contrast  observa¬ 
tions.  Transmittance  depends  on  the  integrated  concentration  along  the  line  of  sight. 

Bleiweiss  and  his  colleagues  at  White  Sands  Missile  Range  supplied  100  ATLAS  transmit¬ 
tance  images  collected  at  0.1-s  intervals  over  a  10-s  period.  As  noted  by  Bleiweiss  et  al.  (1991) 
many  assumptions  are  required.  They  assume  single  scattering  and  a  spatially  uniform  medium 
so  that  transmission  through  the  smoke  Tc  can  be  expressed  in  terms  of  source  and  cloud  radi¬ 
ance  Ls  and  Lc,  and  received  radiance  Lr  as  follows: 

Tc  =  (Lr-Lc)/(Ls-Lc)  (50) 
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FIGURE  6 


AEROSOL  BACKSCATTER  CROSS  SECTIONS  THROUGH  A  PLUME 
APPROXIMATELY  8  KM  DOWNWIND  OF  THE  KINCAID,  ILLINOIS  POWER— 
JULY  1980,  APPROXIMATELY  0950  CST 


(a)  LIDAR  CROSS  SECTION  117 
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FIGURE  6  AEROSOL  BACKSCATTER  CROSS  SECTIONS  THROUGH  A  PLUME 

APPROXIMATELY  8  KM  DOWNWIND  OF  THE  KINCAID,  ILLINOIS  POWER- 
JULY  1980,  APPROXIMATELY  0950  CST  (continued) 


(b)  UDAR  CROSS  SECTION  118 
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FIGURE  6  AEROSOL  BACKSCATTER  CROSS  SECTIONS  THROUGH  A  PLUME 

APPROXIMATELY  8  KM  DOWNWIND  OF  THE  KINCAID,  ILLINOIS  POWER- 
JULY  1980.  APPROXIMATELY  0950  CST  (concluded) 


(c)  LIDAR  CROSS  SECTION  119 
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Blciwciss  et  al.  (1991)  assume  that  cloud  transmittances  Tc  along  nearby  lines  of  sight  with 
contrasting  background  radiance  are  equal  and  that  there  are  no  radiation  sources  off-line  so  that 
Lc  can  be  estimated  from  the  source  and  received  radiances.  If  the  two  lines  of  sight  and  their 
corresponding  background  and  received  radiances  arc  designated  by  subscripts  1  and  2,  Eq.  (50) 
can  be  rearranged  to  give 


L  L  -  L  L 
SI  R2  S2  R1 


(L  -  L  )  (L  -  L  ) 
^  SI  S2'^  R2  Rr 


(51) 


Data  presented  by  Bleiweiss  et  al.  (1991)  show  that  the  ratio  in  Eq.  (51)  remains  nearly  con¬ 
stant  during  an  experiment,  thereby  providing  a  means  for  estimating  L^,  and  from  that  the 
transmittance  distribution.  However,  the  data  exhibit  a  strong  peak  in  the  horizontal  wave  num¬ 
ber  that  corresponds  to  about  three  pixels  (approximately  1  m).  It  is  not  clear  whether  this  is  an 
artifact  of  the  data- reduction  technique,  but  data  at  small  spatial  resolutions  may  not  provide  reli¬ 
able  results  because  of  it. 

Figure  7  shows  three  128  x  128  pixel  subsections  of  the  transmittance  images  in  the 
sequence  supplied  by  Bleiweiss.  Here,  bright  areas  indicate  higher  transmittance,  or  less  aerosol 
along  the  path.  They  correspond  to  times  i.5, 5.5,  and  9.5  s  from  the  beginning  of  the  sequence. 
According  to  information  supplied  by  Bleiweiss  with  the  data,  the  “smoke”  was  an  aluminum 
aerosol  released  from  three  generators  about  500  m  upwind  from  the  part  of  the  plume  shown  in 
Figure  7.  The  instrument  was  about  500  m  from  the  plume.  The  measurements  were  made  at 
1 831  EST,  16  May  1990  at  Eglin  Air  Force  Base,  Florida.  The  winds  at  10  m  were  about  4.4 
ms“l  from  the  south.  The  atmosphere  was  neutral  to  slightly  stable. 


3.  Vector  Data 

A  unique  set  of  atmospheric  wind  observations  was  made  available  to  us  by  Schneider 
(1991).  The  data-measurement  and  reduction  techniques  are  described  in  detail  in  the  above 
cited  reference.  Very  briefly  these  data  were  collected  by  observing  the  motions  of  small  alu¬ 
minum  dipoles  (chaff)  with  two  Doppler  radar  systems  separated  by  about  16  km.  The  two 
radars  provided  data  from  a  volume  about  9  x  9  km  in  the  horizontal  and  2  km  in  the  vertical 
direction  at  a  rate  of  about  one  complete  volume  measured  every  2  min.  The  initial  radial  veloc¬ 
ity  measurements  were  interpolated  to  a  Cartesian  grid  with  200-m  spacing  in  all  three  directions. 
The  grid  is  oriented  with  the  Y  direction  toward  north.  The  radial  velocity  data  were  smoothed 
to  remove  4(K)-m  fluctuations.  The  data  were  also  linearly  interpolated  to  a  common  time 
between  two  successive  volume  scans.  The  u  and  v  components  are  extracted  from  the  two 
radial  velocities  at  each  grid  point,  first  by  ignoring  the  w  component.  The  w  component  is  then 
estimated  from  integration  of  the  continuity  equation  and  used  to  correct  the  first  u  and  v  esti¬ 
mates.  Schneider  (1991)  found  that  the  integration  of  the  continuity  equation  gives  more  reliable 
estimates  of  w  when  performed  in  a  coplanar,  cylindrical  coordinate  system  oriented  along  the 
radar  baseline. 
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FIGURE  7 


A  SEQUENCE  OF  INFRARED  TRANSMITTANCE  IMAGES  THROUGH  AN  ALUMINUM 
AEROSOL  PLUME  (continued) 

Courtesy  Bleiweiss  et  al,.  1991 
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FIGURE  7  A  SEQUENCE  OF  INFRARED  TRANSMITTANCE  IMAGES  THROUGH  AN  ALUMINUM 
AEROSOL  PLUME  (concluded) 

Courtesy  Bleiweiss  et  al.,  1991 


The  data  provided  by  Schneider  were  in  the  form  of  the  u,  v,  and  w  components  at  the  grid 
points.  It  should  be  noted  that  the  chaff  tracked  by  the  radar  was  not  uniformly  distributed  so 
that  there  are  some  “holes”  in  the  data.  Furthermore,  because  the  volumes  scanned  by  the  radars 
were  such  that  there  was  no  overlap  in  some  comers,  those  grid  points  are  also  without  data. 
Nevertheless,  these  data  are  complete  enough  that  I  began  to  consider  applications  to  three- 
dimensional  vector  fields  of  the  fractal  dimension  estimation  methods  discussed  earlier.  As  of 
the  writing  of  this  report,  we  have  not  fully  developed  the  necessary  methodologies.  The  dis¬ 
cussion  of  the  dual  Doppler  wind  data  is  included  to  demonstrate  that  there  are  data  to  which  the 
methods  described  later  can  profitably  be  applied. 

C.  APPLICATIONS  OF  STANDARD  FRACTAL  ANALYSIS  TECHNIQUES 

1.  Fourier-Synthesized  Fractal  Arrays 

The  data  shown  in  Figure  4  were  generated  to  demonstrate  the  degree  to  which  the  different 
fractal  dimension  estimation  techniques  can  recover  the  “true”  fractal  dimension.  In  this  case, 
the  data  are  known  to  conform  to  the  definition  of  a  fractal,  but  (as  we  shall  see)  data  collected  in 
the  atmosphere  do  not  behave  ideally.  We  would  expect  that,  inasmuch  as  the  test  data  were 
generated  by  Fourier  synthesis,  the  Fourier  methods  for  estimating  fractal  dimension  would  work 
well.  This  proves  to  be  the  case  as  can  be  seen  in  Figure  8.  The  three  panels  in  Figure  8  show 
the  power  sjjectra  (averaged  according  to  radial  distance  from  the  origin  in  the  Fourier  plane). 

As  expected,  the  slopes  of  the  best-fit  lines  correspond  to  the  correct  fractal  dimension  within  2% 
in  the  worst  case.  This  is  indicative  of  the  degree  to  which  fractal  dimension  can  be  recovered 
under  the  best  of  circumstances  for  an  array  of  this  size  (256^). 

Figure  9  shows  the  results  obtained  when  the  box -counting  approach  was  applied  to  these 
same  data.  In  this  case,  the  data  were  treated  as  a  three-dimensional  array.  A  240  x  240  subsec¬ 
tion  was  taken  from  the  larger  array.  The  values  at  each  grid  point  were  scaled  so  that  the  largest 
values  were  also  5  240.  In  this  case,  lines  with  slopes  corresponding  to  the  exact  fractal  dimen¬ 
sions  have  been  drawn.  They  show  good  agreement  with  the  points  over  most  of  the  range,  but 
the  points  corresponding  to  the  smaller  sizes  tend  to  show  too  few  points.  If  best-fit  lines  are 
calculated  for  all  but  several  of  the  upper  points,  the  slopes  agree  within  a  Tew  percent  of  the 
exjjected  values. 

Figure  10  shows  the  results  obtained  when  the  multiresolution  figure  analysis  is  applied 
using  the  four-orientation  edge  detector  discussed  earlier  [see  Eq  (13)].  Without  any  formal 
quantification  of  the  congruence  of  the  curves  for  different  thresholds,  this  technique  does  not 
have  the  sensitivity  of  the  others.  It  should  not  be  too  difficult  to  develop  a  method  for  measur¬ 
ing  the  overall  discrepancies  among  the  curves,  but  that  has  not  yet  been  done.  Nevertheless,  it  is 
apparent  that  the  center  panel  of  each  figure,  which  corresponds  to  the  correct  fractal  dimension, 
has  the  most  nearly  coincident  curves. 
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FIGURE  8  POWER  SPECTRA  FOR  SYNTHESIZED  BROWNIAN  FRACTALS 
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FIGURE  9  BOX-COUNTING  ANALYSIS  OF  SYNTHESIZED 


BROWNIAN  FRACTALS 
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2.  Lidar  Data 


In  principle,  there  is  no  reason  why  the  same  analysis  techniques  used  in  the  preceding  sec¬ 
tion  cannot  be  applied  to  any  two-dimensional  array  of  numbers,  such  as  those  representing  the 
lidar  images  of  Figure  6.  The  presumption  is  that  if  the  array  has  the  properties  associated  with 
fractals,  then  that  can  be  determined  along  with  the  fractal  dimension.  The  presence  of  noise  and 
errors  in  the  data  makes  the  problem  more  difficult.  Figure  1 1  shows  the  power  sp)ectra  obtained 
from  the  lidar  images  in  Figure  6.  Figure  1 1  also  shows  the  best-fit  straight  lines  and  their 
slopes.  The  fractal  dimensions  inferred  from  these  slopes  are  generally  >  3,  which  seems  unreal¬ 
istic.  The  box-counting  estimates  shown  in  Figure  12  suggest  a  value  between  2.5  and  2.6,  at 
least  over  the  range  of  larger  sizes  where  the  straight  lines  fits  are  appropriate.  Dashed  straight 
lines  of  slope  -2.9  (corresponding  to  a  fractal  dimension  of  2.55)  have  been  entered  in  Figure  1 1. 
It  is  apparent  that  they  are  in  reasonable  agreement  with  the  box-counting  results  for  the  lower 
wave  numbers  (larger  wavelengths).  It  is  not  certain  whether  the  aerosol  backscatter  is  not  scal¬ 
ing  over  the  entire  range  of  sizes,  or  whether  there  is  random  noise  in  the  data  on  a  smaller  scale. 
Data  uncorrected  for  attentuation  would  be  expected  to  have  somewhat  greater  power  at  wave 
numbers  associated  with  large-scale  features  (Uthe,  1991  personal  communication).  More  recent 
data  from  improved  lidar  systems  and  correction  for  attenuation  may  resolve  this  question. 

Figure  13(a)  is  an  example  of  the  application  of  the  multidimensional  feature- analysis  tech¬ 
nique  to  the  lidar  plume  imagery.  This  figure  is  based  on  the  occurrence  of  “edges”  in  the 
imagery.  For  purposes  of  comparison.  Figure  13(b)  shows  a  similar  analysis  of  the  same  image, 
based  on  the  occurrence  of  peaks  [Eq.  (14)].  In  both  cases,  it  is  obvious  that  the  fine  resolution 
(one  cell,  or  pixel)  counts  do  not  scale  the  same  as  the  coarser  features.  This  is  certainly  consis¬ 
tent  with  the  findings  from  the  spectral  and  box-counting  analyses.  It  also  appears  that  the  two 
features — edges  and  boxes — have  different  scaling  properties.  Figure  13  shows  the  results  for 
the  three  values  of  the  exponent  H  that  were  judged  to  give  the  nearest  congruence  of  the  curves. 
Figure  13(a)  suggests  that  the  best  estimate  of  fractal  dimension  derived  from  the  edge  analysis 
would  be  on  the  order  of  2.6  or  2.7,  reasonably  consistent  with  the  other  methods.  The  peak 
analysis  in  Figure  13(b)  gives  an  estimate  nearer  2.9,  considerably  different  from  the  other  meth¬ 
ods.  Although  the  congruence  of  the  graphs  in  the  figures  is  poor,  especially  compared  to  the 
ideal  cases  shown  in  Figure  10,  the  fact  that  the  two  features  give  very  different  results  suggests 
that  the  scaling  properties  may  depend  on  the  nature  of  the  feature  selected.  This  in  turn  suggests 
that  appropriate  features  may  not  have  been  selected.  This  will  be  discussed  further  in  a  later 
section  of  this  report. 

3.  Transmittance  Imagery 

Applications  of  the  fractal  analysis  techniques  were  least  successful  for  the  transmittance 
imagery.  This  may  be  a  result  of  the  data-reduction  techniques  that  were  used.  As  noted  earlier, 
there  appears  to  be  an  artifact  in  the  data  that  produces  peaks  in  the  horizontal  sp)ectrum  at  wave¬ 
lengths  on  the  order  of  a  meter.  Bleiweiss  et  al.  (1991)  discuss  other  data-reduction  approaches 
in  their  paper  that  do  not  rely  on  the  same  assumptions.  Future  studies  should  compare  results 
obtained  from  images  using  different  data  reduction  methodologies. 
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FIGURE  1 1  POWER  SPECTRA  FOR  LIDAR  PLUME  IMAGES 
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FIGURE  13  MULTIRESOLUTION  ANALYSIS  OF  FIRST  LIDAR  PLUME  IMAGE  (concluded) 
(b)  PEAK  ANALYSIS 
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In  spite  of  problems  with  the  data,  some  results  were  obtained.  Figure  14  shows  the  power 
spectra  for  the  three  transmittance  images  in  Figure  7.  There  is  considerable  variation  among  the 
images.  It  should  also  be  noted  that  if  the  larger  wave  number  (smaller  wavelength)  data  points 
were  ignored,  the  slope  of  the  lines  would  be  steeper,  giving  lower  fractal  dimensions.  There  is 
reasonably  good  agreement  between  the  estimates  of  fractal  dimension  from  the  spectral  slope 
and  those  from  the  box-counting  method  shown  in  Figure  15  for  the  first  and  last  of  the  images. 
The  other  image  is  estimated  to  liave  a  lower  dimension  from  the  spectral  slope  than  from  the 
box  counting. 

As  with  the  lidar  images,  the  transmittance  images  do  not  give  very  good  results  with  the 
multidimensional  feature-analysis  approach.  The  origin  of  the  problem  is  not  clear.  It  may  be 
the  use  of  inappropriate  features,  or  the  data-rcduction  techniques  mentioned  eailier.  Because 
the  transmittance  data  are  appropriate  for  the  study  of  questions  of  considerable  interest  to  the 
Army,  every  effort  should  be  made  in  the  future  to  use  these  data  to  determine  the  nature  of 
inhomogeneities  in  the  transmittance  through  smoke  plumes. 

D.  EXTENSIONS  OF  MULTIRESOLUTION  FEATURE  ANALYSIS 
1.  Vectors 

a.  Two-Dimensional  Fields 

Scalar  fields  have  been  the  focus  of  the  discussion  to  this  point,  but  multircsolution  feature 
analysis  can  also  be  applied  to  vector  fields.  Two  possible  approaches  exist: 

•  Apply  any  of  the  techniques  to  some  scalar  property  of  the  vector  field,  such  as  diver¬ 
gence  or  the  vertical  component  of  vorticity. 

•  Define  vector  features  and  apply  the  multiresolution  methodology  directly  to  the  vector 
field. 

Some  of  the  conventional  finite-difference  approximations  for  vector  field  properties  can  be  rep¬ 
resented  by  features  like  that  represented  in  Figure  16.  It  shows  the  finite-difference  operator  for 
the  vertical  component  of  vorticity  expressed  as  the  sum  of  two  features  (like  those  presented 
earlier  as  matrices  to  be  superimposed  on  the  data  field  to  define  a  sum  of  products):  One  of 
these  “templates”  is  applied  to  the  westerly  (u)  component,  the  other  to  the  southerly  (v)  compo¬ 
nent.  Ordinary  scalar  feature  detectors  can  then  be  applied  to  the  resulting  feature  intensity  field, 
although  the  interpretation  may  not  be  straightforward.  Figure  17  shows  a  flow  field  that  would 
be  associated  with  an  “edge”  in  the  vorticity  field. 

Scalar  properties  of  the  vector  field  need  not  be  used;  interpretation  of  the  results  may  be 
much  easier  if  vector  features  are  defined  and  applied  directly  to  the  field.  Figure  18  shows  two 
particularly  interesting  vector  features.  The  vector  templates  are  applied  in  the  same  way  as  the 
scalar  templates,  except  that  the  arithmetic  product  is  replaced  by  the  scalar  (dot)  product  of  each 
feature  vector  and  the  corresponding  vector  in  the  field  being  analyzed.  If  the  vortex  feature 
template  in  Figure  18(b)  is  applied  to  variously  smoothed  vector  fields,  then  we  are  actually 
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FIGURE  14  POWER  SPECTRA  FOR  ATLAS  TRANSMISSOMETER  IMAGES 
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U  V 

FIGURE  1 6  FEATURE  BASED  ON  THE  FINITE-DIFFERENCE  APPROXIMATION  FOR 
THE  VERTICAL  COMPONENT  OF  VORTICITY 

pursuing  the  quantitative  meaning  of  Richardson’s  doggerel  that  “Big  whorls  have  little 
whorls. . (Richardson,  1922,  as  quoted  by  Mandelbrot,  1983).  The  whorls  or  vortices  in 
Figure  18(b)  are  two-dimensional  features,  that  might  be  associated  with  the  feature  in  Figure 
18(a)  oriented  in  a  perpendicular  plane. 

b.  Three-Dimensional  Fields 

One  of  the  mechanisms  suggested  for  the  transfer  of  energy  down  the  turbulent  cascade  is 
by  the  stretching  of  vortex  filaments.  Figure  18(a)  shows  a  flow  pattern  that  produces  stretching 
along  the  axis  between  the  upper  left  and  lower  right  comers.  If  this  transfer  mechanism  is 
important,  then  the  feature  shown  in  Figure  18(b)  should  have  large  (positive  or  negative)  values 
on  a  smaller  scale  in  the  plane  normal  to  the  axis  of  stretching.  This  relationship  needs  to  be 
studied  in  applications  to  observed  and  simulated  flow  fields.  As  just  described,  the  method 
would  be  applied  to  two-dimensional  features  in  two  steps.  It  would  be  worthwhile  looking  at 
three-dimensional  features,  but  to  do  so  would  be  difficult.  In  the  example  just  discussed,  the 
features  in  Figure  18(a)  will  be  of  a  smaller  scale  than  those  in  Figure  1 8(b),  and  they  will  be  ori¬ 
ented  normal  to  the  stretching  axis. 

Figure  19  shows  an  example  of  a  three-dimensional  vector  feature  that  corresponds  to  verti¬ 
cal  stretching  of  an  eddy  in  the  horizontal  plane.  In  this  case,  the  value  assigned  for  the  feature  at 
the  center  point  would  be  the  sum  of  the  scalar  products  of  the  feature  vectors  with  the  corre¬ 
sponding  vector  in  the  field  being  analyzed.  Once  the  feature  strengths  have  been  determined  for 
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edge  in  vorticity  field 


(a)  (b) 

FIGURE  18  TWO-DIMENSIONAL  VECTOR  FEATURES  FOR 
TWO  DIFFERENT  TYPES  OF  SINGULARITY 

(a)  HYPERBOUC  POINT 

(b)  COUNTERCLOCKWISE  VORTEX 


FIGURE  17 


AN  EDGE  FEATURE  AND  A  VELOCITY  FIELD  THAT  PRODUCES 
AN  EDGE  IN  THE  VORTICITY  FIELD 


FIGURE  19  EXAMPLE  OF  A  THREE-DIMENSIONAL  VECTOR  FEATURE 

variously  smoothed  vector  fields,  the  information  necessary  for  multiresolution  analysis  is  avail¬ 
able.  The  dual  Doppler  radar  wind  data  discussed  earlier  should  provide  a  good  basis  for  case 
studies.  Although  the  number  of  data  points  in  any  one  field  may  be  too  small,  grouping  the  data 
should  give  sufficiently  large  samples  for  analysis. 


2.  Empirically  Defined  Features 

One  problem  that  arises  in  the  use  of  vector  features  (and  to  a  lesser  extent,  scalar  features 
as  well)  is  that  one  quickly  runs  out  ideas  for  features.  Furthermore,  unless  one  is  careful,  there 
may  be  considerable  redundancy  of  information  among  the  members  of  any  set  of  features  that  is 
arbitrarily  defined  (i.e.,  they  may  be  linear  combinations  of  one  another).  Thus,  to  carry  this 
approach  the  logical  next  step,  it  would  be  desirable  to  have  an  approach  to  the  definition  of  fea¬ 
tures  that  was  based  on  the  characteristics  of  the  data  themselves,  and  that  provided  features  that 
were  independent  of  one  another.  It  would  also  be  desirable  to  have  a  measure  of  the  relative 
importance  of  the  different  features. 

In  the  case  of  scalar  features,  an  approach  that  is  widely  known  is  available.  For  example, 
Lorenz  ( 1 956)  represented  the  patterns  of  variability  of  atmospheric  pressure  at  64  stations  in  the 
United  States  as  a  linear  combination  of  independent  patterns  of  variability,  which  he  referred  to 
as  Empirical  Orthogonal  Functions  (EOF).  Lumley  (1967,  1980)  suggested  that  a  similar  kind 
of  analysis  could  be  used  to  extract  coherent  structures  from  turbulent  flows.  Ludwig  and  Byrd 
(1980)  also  applied  the  concept  to  vector  fields.  They  identified  patterns  of  variability  in  the 


52 


inputs  used  for  a  wind  model  in  order  to  simplify  the  resulting  calculations.  Sirovich  (1988) 
describes  the  analysis  of  turbulent  flows  by  a  similar  procedure.  According  to  Sirovich 
(personal  communication,  1991),  the  patterns  of  variability  can  be  identified  without  as  much 
calculation  as  is  required  by  the  covariance  matrix  diagonalization  approach  described  by  Lorenz 
(1956),  and  adapted  for  use  with  vectors  by  Ludwig  and  Byrd  (1980).  The  newer  approach  has 
not  been  tried  yet. 

The  approach  that  has  been  pursued  here  uses  small  3x3x3  subsections  of  the  observed 
wind  field.  Subtracting  the  center  vector  from  each  of  the  surrounding  26  vectors  in  the  subsec¬ 
tion  gives  a  difference  vector  (Av)  for  each  of  the  points  in  the  subsection,  i.e.. 
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where  UMk,vijk  and  wjjk  are  vector  components  at  point  ijk  (222  is  at  the  center  of  the  subsection) 
in  the  m^  subsection.  The  three-dimensional  array  of  AVs  shows  the  local  pattern  of  motion 
about  the  center  point.  If  we  have  this  array  and  the  vector  V222  at  the  center  point,  then  we  can 
reconstruct  the  26  surrounding  vectors.  Looking  ahead  to  possible  applications,  we  can  envision 
that  we  have  modeled  the  vectors  on  a  coarse  array  (corresponding  to  the  center  point  values). 
Then,  if  we  have  some  way  of  estimating  the  array  of  AVs,  we  can  obtain  the  field  with  finer 
resolution. 

Next,  we  determine  the  deviations  AV'  about  the  means  by  subtracting  the  average  (indi¬ 
cated  by  the  overbar)  over  the  same  relative  points  in  all  N  subsections  from  individual  aV’s, 
i.e.,: 


(^''ijO  =  .  (53) 

m  m  ' 

This  relative  variability  about  the  center  point  is  described  by  (3  x  3  x  3  grid  points)  x  3 
vector  components.  These  81  numbers  are  treated  as  the  components  of  a  column  vector  describ¬ 
ing  the  “state”  of  the  subsection.  A  matrix  of  these  “state  vectors”  is  multiplied  by  its  transpose 
to  give  the  covariance  matrix  for  the  complete  set  of  state  vectors.  The  eigenvectors  of  this 
matrix  that  account  for  the  most  variance  in  the  individual  patterns  can  be  used  as  the  “features” 
in  the  modified  multiresolution  feature  analysis  methodology.  These  eigenvectors  are  the  empir¬ 
ical  orthogonal  functions  (EOFs)  that  will  be  used  for  subsequent  feature  analysis. 

The  same  approach  can  as  easily  be  applied  to  scalar  fields,  differing  only  in  that  there  is  but 
one  number  associated  with  each  grid  point,  so  that  it  is  computationally  feasible  to  look  at  larger 
patterns  of  variability  with  scalars.  Identification  of  the  major  patterns  of  variability  in  scalar 
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fields  and  how  they  scale  with  size  could  probably  be  used  to  guide  smoke  plume  simulation 
methodologies  like  that  recently  described  by  Hoock  (1991).  He  redistributes  smoke  material  to 
smaller  size  elements  based  on  the  space -fillingness  derived  from  an  estimate  of  the  fractal 
dimension.  It  is  easy  to  see  that  the  redistribution  could  be  guided  by  the  observed  patterns  of 
variability  and  their  relative  frequencies,  rather  than  by  purely  random  redistribution. 

FORTRAN  programs  have  been  written  to  determine  the  most  important  features  in  three- 
dimensional  vector  arrays  and  two  dimensional  scalar  arrays.  The  programs  were  applied  to 
determine  the  vector  features  for  the  dual  Doppler  radar  observations  of  winds  near  Boulder, 
Colorado  for  1246  MST  on  22  June  1984  (Schneider,  1991).  The  average  variability  about  the 
center  point  of  3  x  3  x  3  subsections  of  the  wind  field  and  the  three  vector  EOFs  that  explain  the 
most  variance  are  shown  in  Figure  20.  The  strong  shear  in  the  wind  field  is  evident  in  the  aver¬ 
age  pattern  of  variability.  The  first  EOF  indicates  a  tendency  for  small-scale  patterns  (on  the 
order  of  400  m  on  a  side)  to  show  strengthening  and  weakening  of  vortex  patterns  tilted  in  the 
approximate  direction  of  the  shear.  The  second  most  important  pattern  of  local  spatial  fluctua¬ 
tion  (explaining  almost  as  much  variance  as  the  first),  is  a  general  strengthening  or  weakening  of 
the  existing  shear.  The  third  EOF  is  somewhat  more  complex  (and  higher  order  EOFs  tend  to  be 
even  more  so);  its  most  important  feature  is  that  it  accounts  for  local  strengthening  (and  weaken¬ 
ing)  of  the  horizontal  shear  in  approximately  the  north-south  direction.  The  results  of  this  appli¬ 
cation  are  encouraging.  It  is  hoped  that  the  techniques  alLded  to  by  Sirovich  (1988)  can  reduce 
the  required  calculations.  In  any  event,  I  expect  to  continue  development  and  application  of 
these  techniques  to  atmospheric  wind  measurements. 
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AVERAGE  VARIATION  ABOUT  THE  CENTER  POINT  AND  THREE  MOST 
IMPORTANT  PATTERNS  OF  VARIATION  IN  AN  OBSERVED  WIND  FIELD 
(PERCENTAGE  OF  EXPLAINED  VARIANCE  IS  SHOWN) 


FIGURE  20 
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IV  FUTURE  RESEARCH 


A.  SUBGRID  FLOW  SIMULATION 
1.  Deterministic  Approaches 

It  is  quite  possible  that  the  EOFs  will  depend  on  the  prevailing  meteorology.  For  example, 
there  is  a  tantalizing  hint  of  this  in  the  first  EOF  shown  in  Figure  20.  Although  it  may  just  be 
happenstance  that  it  has  an  axis  that  tilts  from  the  vertical  in  the  direction  of  the  large-scale  shear, 
it  would  not  be  surprising  to  find  that  the  EOFs  are  affected  by  such  things  as  the  large-scale  cir¬ 
culation,  thermal  stratification,  and  so  forth.  Furthermore,  the  most  important  EOFs  may  exhibit 
important  changes  from  one  level  of  resolution  to  the  next.  Such  variability  would  complicate 
matters,  but  should  not  prevent  the  kinds  of  analyses  discussed  below  from  being  performed. 

It  should  somehow  be  possible  to  characterize  the  nature  and  strength  of  the  most  important 
features  at  a  small  scale  from  information  available  at  coarser  scales.  If  the  field  is  fractal,  then  it 
should  at  least  be  possible  to  define  the  probability  distribution  of  feature  strengths,  opening  the 
way  to  probabilistic  simulations  of  small-scale  features.  Of  course,  it  will  be  easier  if  some  rea¬ 
sonable  correlations  can  be  found  between  feature  strengths  at  on  scale  and  those  at  smaller  reso¬ 
lutions.  This  would  open  the  door  to  a  simpler,  more  “deterministic”  relationship  between 
scales.  There  are  two  possibilities,  but  in  either  case  the  process  begins  by  defining  the  several 
most  important  EOFs  for  a  combined  set  of  wind  fields  that  correspond  to  similar  meteorological 
conditions.  The  EOFs  will  be  defined  for  each  resolution  and  for  the  combined  data  set  compris¬ 
ing  of  all  the  different  smoothings.  Then,  either  the  EOFs  for  the  combined  resolutions  closely 
re.semble  those  derived  for  each  of  the  smoothings,  or  they  differ  substantially  from  one  another. 
In  the  first  case,  we  choose  one  set  of  EOFs  from  those  that  have  been  calculated.  Initially,  this 
will  probably  be  a  subjective  choice,  but  an  objective  selection  scheme  based  on  the  inner  prod¬ 
ucts  of  the  EOFs  (which  are  actually  large  vectors)  should  not  be  difficult  to  devise.  The  next 
step  is  to  use  the  most  important  EOFs  as  features  and  determine  their  strength  at  the  grid  pioints 
for  each  smoothing  in  each  of  the  wind  fields  of  the  set.  These  feature  strengths  can  then  be  used 
directly  in  the  scheme  adapted  from  Jones  et  al.  (1991)  to  estimate  fractal  dimensions  and 
determine  whether  different  values  are  associated  with  different  features.  These  same  data  also 
provide  the  information  necessary  for  developing  deterministic  or  probabilistic  relationships  as 
described  later. 

If  the  EOFs  are  substantially  different  for  different  resolutions,  then  we  can  still  use  the 
EOFs  derived  from  the  combined  smoothings  as  a  basis  for  estimating  fractal  dimensions.  How¬ 
ever,  for  estimating  smaller-scale  circulations,  it  will  be  better  to  calculate  feature  strengths  for 
each  smoothed  field  based  on  the  EOFs  that  have  been  found  to  explain  the  most  variance  for 
that  particular  resolution.  Again,  this  will  lead  to  a  set  of  feature  values  associated  with  each  grid 
point  in  each  of  the  smoothed  fields.  However,  those  feature  strengths  will  refer  to  different 
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features  at  the  different  resolutions.  This  will  not  necessarily  present  a  problem  when  we  search 
for  deterministic  or  probabilistic  relationships  between  features  at  different  scales,  but  it  is  likely 
to  introduce  some  difficulties  in  practical  applications. 

One  approach  to  the  analysis  would  begin  by  calculating  the  linear  correlations  between  the 
feature  strength  (for  each  of  the  most  important  small-scale  features)  at  points  on  the  smallest 
scale  grid,  and  the  feature  values  for  the  most  important  features  at  nearby  points  on  the  next 
coarser  grid.  There  are  some  limitations  on  the  choice  of  variables  that  are  imposed  by  the 
geometry,  i.e.,  the  coarse  grid  points  for  which  feature  values  can  be  calculated  will  all  be  farther 
removed  from  the  edges  of  the  domain  than  the  finer  scale  points.  Nevertheless ,  it  is  likely  that 
we  will  end  up  with  a  correlation  matrix  of  about  4  (important  fine-scale  features)  x  108 
(representing  four  features  for  each  of  the  nearest  27  points).  If  the  linear  correlation  table  does 
not  have  any  high  correlations,  it  would  probably  be  worthwhile  to  repeat  the  process  using  rank 
correlations  in  case  there  are  significant  monotonic,  but  nonlinear,  relationships  in  the  data. 

If  significant  correlations  are  found,  their  distribution  through  the  table  may  be  very  enlight¬ 
ening.  For  example,  if  the  correlations  are  higher  for  feature  values  at  more  distant  points  than 
for  the  values  at  the  central  point,  it  suggests  that  the  transfer  of  motion  to  smaller  scales  involves 
small-scale  features  “attached”  to  the  larger  features,  rather  than  “embedded”  in  them.  Correla¬ 
tions  between  different  kinds  of  features  could  provide  hints  about  the  physical  processes 
involved.  For  example  correlation  between  small  scale  “eddy-like”  features  and  larger-scale 
“stretching”  features  might  be  evident  if  the  oft-suggested  cascade  via  vortex  stretching  is 
important. 

If  some  instances  of  high  correlations  are  found,  then  those  coarse  features  that  are  most 
highly  correlated  with  the  smaller-scale  feature  values  would  be  used  as  input  to  an  optimized 
multiparameter  regression  scheme  to  provide  a  functional  relationship  between  the  values  of 
coarse-scale  features  and  those  of  the  most  important  smaller-scale  ones.  If  the  relationships 
were  generally  valid  for  all  cases  with  similar  meteorological  conditions,  then  the  regression 
results  would  provide  all  the  required  information  for  an  extrapolation  to  smaller  scales.  The 
larger-scale  feature  intensities  can  be  derived  directly  from  the  large-scale  grid  values,  whether 
obtained  by  simulation  or  observation.  Those  feature  values  in  turn  are  used  to  estimate  the  fea¬ 
ture  values  for  the  next  smaller  scale,  which  serves  as  basis  for  estimating  small-scale  values 
from  the  center  point  value  (which  is  at  a  coarse  grid  point),  average  variations  about  that  point, 
and  a  linear  combination  of  the  features,  using  their  estimated  strengths  as  the  coefficients. 

2.  Probabilistic  Approaches 

If  there  are  no  strong  correlations  between  feature  strengths  at  the  smallest  scale  and  those 
at  the  next  largest,  the  problem  is  much  more  difficult,  but  there  are  at  least  two  possible 
approaches  that  can  be  tried.  Both  are  related  to  the  probability  of  occurrence  of  feature 
strengths  (or  combinations  of  feature  strengths)  on  the  small  scale  to  those  on  the  larger  scale. 

The  fact  that  little  or  no  correlation  was  found,  combined  with  the  presumed  fractal  nature  of  the 
fields,  suggests  the  use  the  scaling  properties  to  define  the  small-scale  feature  intensity  probabil¬ 
ity  distributions  and  use  them  in  a  Monte  Carlo  scheme  to  define  a  small-scale  feature  strength 
for  each  important  feature  around  each  coarse  grid  point.  The  orthogonality  of  the  features 
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guarantees  that  their  strengths  are  uncorrelated,  so  the  probability  distributions  can  be  used 
independently  without  worrying  about  possible  cross  correlations.  Once  a  feature  strength  has 
been  defined  for  each  feature  type  at  each  grid  point,  the  process  of  constructing  the  smaller  scale 
values  proceeds  just  as  described  above  for  the  deterministic  case. 

If  no  strong  correlations  arc  found  between  the  strengths  of  individual  feature  types  at  the 
different  scales,  it  would  be  worth  exploring  the  possibility  of  using  “clusters”  of  feature 
strengths  to  stratify  the  analysis.  One  of  the  standard  cluster  analysis  programs  would  be  used  to 
identify  groupings  of  feature  strengths  at  the  coarse  and  fine  scales.  Each  coarse  grid  point 
would  be  identified  with  one  of  the  coarse  clusters  and  with  one  of  the  fine  clusters.  Then,  con¬ 
tingency  tables  would  be  constructed  to  show  the  probability  of  occurrence  for  each  fine-scale 
cluster,  given  the  occurrence  of  a  specified  coarse  cluster.  These  contingency  tables  could  be 
used  in  a  Monte  Carlo  scheme  to  select  a  fine-scale  cluster  for  each  grid  point,  given  the  coarse 
cluster  that  was  appropriate  to  that  point.  We  could  calculate  an  average  value  for  each  feature 
strength,  based  on  the  members  of  the  cluster,  and  use  those  feature  strengths  as  described  before 
to  construct  a  local  fine  scale  field.  In  essence,  if  this  approach  were  applied  to  scalars,  it  would 
be  a  refinement  of  that  described  by  Hoock  (1991).  He  selects  the  smaller  cells  in  which  material 
is  to  appear  randomly  (with  the  average  number  of  such  cells  defined  by  the  fractal  dimension), 
rather  than  by  having  preferred  patterns  of  distribution. 

With  incomplete  approximation  (not  all  EOFs  used),  it  is  likely  that  the  local  vector  (or 
other  fields  will  violate  various  laws  of  physics — especially  if  full  set  of  dynamic  variables  are 
included.  Thus,  it  may  be  necessary  to  adjust  the  first-guess  field  to  be  compatible  with  govern¬ 
ing  equations.  Possibilities  include  iterative  relaxation  schemes  and  variational  calculus 
approach  to  get  minimum  adjustment. 

B.  fractal/scaling  properties 

The  multiresolution  feature-analysis  methodology  requires  the  calculations  of  “feature 
strength”  fields  at  all  scales.  These  are  scalar  fields  (regardless  of  whether  the  original  field  was 
scalar  or  vector)  that  can  be  used  as  input  to  any  of  the  other  fractal  dimension  calculation 
schemes,  which  would  provide  a  link  to  other  studies.  It  would  also  serve  as  another  method  for 
examining  whether  different  aspects  of  the  field  are  distributed  in  space  with  different  scaling 
properties.  This  would  not  be  surprising,  as  the  different  patterns  represent  different  aspects  of 
the  flow  (or  scalar  distribution).  As  an  example,  it  was  shown  earlier  that  the  vertical  component 
of  vorticity  can  be  expressed  as  a  vector  feature.  Its  distribution  might  well  be  expected  to  be 
different  from  that  of  other  scalars,  such  as  speed,  divergence,  or  other  flow  characteristics. 

One  interesting  question  that  may  need  to  be  addressed  is  the  following:  Does  the  variance 
explained  by  the  EOFs  come  from  many  common  “events”  of  modest  size,  or  from  a  few 
extreme  cases?  The  former  would  be  more  useful  for  ordinary  fluid  flow  parameterization  appli¬ 
cations.  The  latter  would  be  more  useful  for  addressing  extreme  value  concerns,  but  it  seems 
unlikely  that  it  would  be  possible  to  acquire  large  enough  data  sets  with  the  required  wide  range 
of  scales. 
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